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Preface 


This is the first annual report to the U.S. Army Medical Research and Material 
Command for the three year project “Advanced Soft Tissue Modeling for Telemedicine 
and Surgical Simulation” supported by grant No. DAMD17-01-1-0673 to The Cleve- 
land Clinic Foundation, to which the NASA Glenn Research Center is a subcontractor 
through Space Act Agreement SAA 3-445. 

The objective of this report is to extend popular one-dimensional (ID) fractional- 
order viscoelastic (FOV) material models into their three-dimensional (3D) equiva- 
lents for finitely deforming continua, and to provide numerical algorithms for their 
solution. The present report is organized into seven chapters and three appendices. 

The first chapter serves as an introduction to the fractional calculus. Algorithms 
for computing fractional derivatives, fractional integrals, fractional-order differential 
equations (FDE’s), and the Mittag-Leffler function (which apprears in analytic solu- 
tions of FDE’s) are provided. 

One of the oldest applications of the fractional calculus is viscoelasticity. Chapter 
two presents an overview of ID FOV. Definitions for the standard FOV fluid and 
the standard FOV solid are put forth along with formulae that are useful in their 
characterization, assuming infinitesimal strains and rotations. 

The third chapter provides an overview of continuum mechanics using body (i.e. , 
convected) tensor fields. Three strain fields are introduced that are measures of 
strain based on changes in: length of line, separation of non-intersecting surfaces, 
and volume of mass. Introduced here for the first time are fractal rates of arbitrary 
tensor fields. Body fields are useful when deriving contitutive equations. 

In the fourth chapter, the body fields defined in the previous chapter are mapped 
into objective, Cartesian, space fields. A useful by-product of field transfer is that 
those spatial fields created by field transfer are frame invariant. Spatial fields are 
useful when solving boundary- value problems. 

The fifth chapter derives isotropic and transverse-isotropic theories for elastic and 
viscoelastic materials by applying a work potential to an integrity basis. Both com- 
pressible and incompressible materials are considered. These theories are derived 
in the body and then transferred into Cartesian space in both the Eulerian and La- 
grangian frames. The tangent modulus is derived for the general theoretical structures 
of elastic and viscoelastic solids. 

A suite of homogeneous experiments used to characterize material models is pre- 
sented in the sixth chapter. The suite includes the homogeneous deformations of: 
shear-free extension (e.g., uniaxial elongation, biaxial extension, pure shear, and di- 
lational compression) and simple shear. The deformation, stress and strain fields 

xiii 
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defined in the prior chapter, along with their various rates, are all quantified for this 
suite of experiments. 

Chapters seven through nine provide elastic and viscoelastic constitutive models 
appropriate for 3D analysis. Chapter seven provides material models for bulk re- 
sponse. Chapter eight will introduce material models for isotropic elastomers, while 
chapter nine will introduce material models for soft biological tissues, which are gen- 
erally transverse isotropic; they will be completed for the second annual report. Both 
classical and fractional-order viscoelastic models are presented. Included are solutions 
for the characterization experiments of chapter six. 

There are three appendices. The first appendix tabulates Caputo fractional deriva- 
tives for a few of the more common mathematical functions. The second appendix 
outlines an automatic procedure for numerical integration that is required by the algo- 
rithm which computes the Mittag-Leffler function. And the third appendix provides 
an efficient scheme for approximating a specific form of the Mittag-Leffler function 
that arises in FOV. 
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Chapter 1 


Fractional Calculus: 
numerical methods 

1.1 Riemann-Liouville Fractional Integral 

In the classical calculus of Newton and Leibniz, Cauchy reduced the calculation of an 
n-fold integration of the function y(x) into a single convolution integral possessing 
an Abel (power law) kernel, 

rx fXn-1 

J n y{x) := / / • • • / y(x 0 ) dx 0 ... dx„_ 2 cte„_ x 

JO Jo Jo ^ 

= <zhjiL ^=w r - y(z ' )dx '' n€N ’ 

where J n is the n-fold integral operator with J° y(x) — y(x), N is the set of positive 
integers, and R + is the set of positive reals. Liouville and Rieraann* analytically 
continued Cauchy’s result by replacing the discrete factorial ( n — 1)! with Euler’s 
continuous gamma function F(n), noting that (n — 1)! = F(n), thereby producing [67, 
15, Eqn. A] 

J °y&- = T jsjl a ' ieR+ - (1 - 2) 

where J a is the Riemann-Liouville integral operator of order a, which commutes (i.e., 
J a J i9 y(x) — jPj a y(x) = J a+ Py(x) V a, 0 G M+). Equation (1.2) is the cornerstone of 
the fractional calculus, although it may vary in its assignment of limits of integration. 
In this report we take the lower limit to be zero and the upper limit to be some 
positive finite real. Actually, a can be complex [102], but for our purposes we only 
need it to be real. 

A brief history of the development of fractional calculus can be found in Ross 
[100] and Miller and Ross [78, Chp. 1]. A survey of many emerging applications of 
the fractional calculus in areas of science and engineering can be found in the recent 
text by Podlubny [86, Chp. 10]. 

’Riemann’s pioneering work in the field of fractional calculus was done during his student years, 
but published posthumous — forty-four years after Liouville first published in the field [100]. 
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1.2 Cap uto- Type Fractional Derivative 

From this single definition for fractional integration one can construct several def- 
initions for fractional differentiation (cf. e.g., [86, 102]). The special operator 
that we choose to use, which requires the dependent variable y to be continuous and 
[a] -times differentiable in the independent variable x , is defined by 

D: y(x) := J M ~ a D M y(x), (1.3) 

such that 

lim D“y(x) = D n y(x) for n E N, (1.4) 

with D\ y(x) = y(x), where [a] is the ceiling function giving the smallest integer 
greater than (or equal to) a, and where a — > n_ means a goes to n from below. 
The operator D n , n € N, is the classical differential operator. It is accepted practice 
to call T>“ the Caputo differential operator of order a, after Caputo [12] who was 
the amoung the first to use this operator in applications and to study some of its 
properties.^ Appendix A presents a table of Caputo derivatives for some of the more 
common mathematical functions. 

The Caputo differential operator is a linear operator 

D*{y + z)(z) = D*y(x) + D“z(x) (1.5a) 

that commutes 

0!Dfy(x) = Dfo; y« = D? +e y(x) V a, 0 € R+ (1.5b) 

if y(x) is sufficiently smooth, and it possesses the desirable property that 

D°c — 0 for any constant c. (1.5c) 

The more common Riemann-Liouville fractional derivative D a , although linear, need 
not commute [86, pg. 74]; furthermore, D a c = I)WjN“ a c = cx~ a /T{l — a), which 
is a function of x\ Ross [100] attributes this startling fact as the main reason why 
the fractional calculus has historically had a difficult time being embraced by the 
mathematics and physics communities. 

* Actually, Liouville introduced the operator in his historic first paper on the topic [67, f 6, Eqn. B]. 
Still, nothing in Liouville’s works suggests that he ever saw any difference between D° = jr a l~ a £)r“l 
and D a = D a ji Q l _a , D a being his accepted definition [67, first formula on pg. 10] — the Riemann- 
Liouville differential operator of order a. Liouville freely interchanged the order of integration and 
differentiation, because the class of problems that he was interested in happened to be a class where 
such an interchange is legal, and he made only a few terse remarks about the general requirements 
on the class of functions for which his fractional calculus works [74], The accepted naming of the 
operator after Caputo therefore seems warrented. 

Rabotnov [90, pg. 129] introduced this same differential operator into the Russian viscoelastic 
literature a year before Caputo’s paper was published. Regardless of this fact, operator £>“ is 
commonly named after Caputo in the current liturature. 


NASA/TM— 2002-21 1914 


2 



The Riemann-Liouville integral operator J a and the Caputo differential operator 
£)“ are inverse operators in the sense that 



a . 

ifc=0 


a e R+, 


( 1 . 6 ) 


with := D k y(0 + ), where [aj is the floor function giving the largest integer less 
than a. The classic ra-fold integral and differential operators of integer order satisfy 
like formulae, viz.: D n J n y(x) = y(x) and J n D n y(x) — y(x) — YllZo fryo+i n € N. 

A word of caution. Fractional derivatives do not satisfy the Leibniz product rule 
of classical calculus. For example, whenever the Caputo derivative is restricted so 
that 0 < a < 1, the Leibniz product rule is given by 


D °{ y x z)(x) 


y(0 + ) z(x) — z(0 + ) , \. . , . 

r(1 r t) x — ^ — + ( D *y)( x ) x < x ) 
+ E (*) ( Jfc " ay K*) x 


(1.7) 


where, unlike the Leibniz product rule for integer-order derivatives, the binomial 
coefficients (“) = (with («) = i ( a e 14 an d k e N) do not 

become zero whenever k > a because a N (i.e., the binomial sum is now of infinite 
extent). A similar infinite sum exists for the Leibniz product rule of the Riemann- 
Liouville fractional derivative (cf. Podlubny [86, pp. 91-97]). 


1.2.1 Integral Expressions 

The Caputo derivative (1.3) can be expressed in more explicit notation as the integral 


D a Mx) = 


r(M 


— r 

- a) J 0 


(x — a:') 0- !-"! 


(D M y)(z') dx ' , 


a, x e 


(1.8a) 


where the weak singularity caused by the Abel kernel of the integral operator is readily 
observed. This singularity can be removed through an integration by parts 


W = r(1 + A- _ a) + £"(*- &') . 

(l-8b) 

provided that the dependent variable y is continuous and (1+ [a]) -times differentiable 
in the independent variable x over the interval of differentiation (integration) [0, x]. In 
(1.8b) the power-law kernel is bounded over the entire interval of integration; whereas, 
in (1.8a) the kernel is singular at the upper limit of integration. 

The two representations of (1.8a) and (1.8b) are quite useful for pen-and-paper 
calculations, but in order to obtain a numerical scheme for the approximation of 
such fractional derivatives, we found it even more helpful to look at yet another 
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representation that seems to have been introduced into this context by Elliott [30]; 
namely, 

B *“ y(x) = f rby/ (x- lr +i y(x ' )dx '' (L8c) 

This representation can also be obtained from (1.8a) using the method of integration 
by parts, but with the roles of the two factors interchanged. The advantage here 
is that the function y itself appears in the integrand instead of its derivative. The 
disadvantage is that the singularity of the kernel is now strong rather than weak, and 
thus we have to interpret this integral as a Hadamard-type finite-part integral. This 
is cumbersome in pen-and-paper calculations but, as we shall see below, it is not a 
problem to devise an algorithm that makes the computer do this job. We provide a 
brief description of such an algorithm in the following pages. For more details, the 
interested reader is referred to [20, 30] and the references cited therein. 

1.3 Caputo-Type FDE’s 

Fractional material models, the subject of this report, are systems of fractional-order 
differential equations (FDE’s) that need to be solved in accordance with appropriate 
initial and boundary conditions. A FDE of the Caputo type has the form 

D*y{x) - f(x,y(x)), a,xeE+, (19a) 

satisfying the (possibly inhomogeneous) initial conditions 

y£+ = D k y( 0 + ), A; = 0,1,..., [aj, (l-9b) 

and whose solution is sought over an interval [0,X], say, where X G JR+. It turns 
out that under some very weak conditions placed on the function f of the right-hand 
side, a unique solution to (1.9) does exist [21]. 

A typical feature of differential equations (both classical and fractional) is the 
need to specify additional conditions in order to produce a unique solution. For the 
case of Caputo FDE’s, these additional conditions are just the static initial condi- 
tions listed in (1.9b), which are akin to those of classical ODE’s, and are therefore 
familiar to us. In contrast, for Riemann-Liouville FDE’s, these additional conditions 
constitute certain fractional derivatives (and/or integrals) of the unknown solution 
at the initial point x = 0 [57], which are functions of x\ These initial conditions are 
not physical; furthermore, it is not clear how such quantities are to be measured from 
experiment, say, so that they can be appropriately assigned in an analysis. I If for no 

*We explicitly note, however, the very recent paper of Podlubny [87] who attempts to give 
highly interesting geometrical and physical interpretations for fractional derivatives of both the 
Riemann-Liouville and Caputo types. These interpretations are deeply related to the questions: 
What precisely is time? Is it absolute or not? And can it be measured correctly and accurately, and 
if so, how? Thus, we are still a long way from a full understanding of the geometric and physical 
nature of a fractional derivative, let alone from an idea of how we can measure it in an experiment, 
but our mental picture of what fractional derivatives and integrals ‘look like’ continues to improve. 
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other reason, the need to solve FDE’s is justification enough for choosing Caputo’s 
definition (i.e., Z?" = jW _a Df Q 1) for fractional differentiation over the more com- 
monly used (at least in mathematical analysis) definition of Liouville and Riemann 
(viz., D a = 


1.4 Numerical Approximations 

1.4.1 Caputo-type Fractional Derivatives 

Unlike ordinary derivatives, which are point functionals, fractional derivatives are 
hereditary functionals possessing a total memory of past states. A numerical algo- 
rithm for computing Caputo derivatives has been derived by Diethelm [20p and is 
listed in Alg. 1.1. Validity of its Richardson extrapolation scheme for 1 < a < 2, 
or one similar to it, has to date not been proven, or disproven. Here y n denotes 
y(rc n ), while yjv represents y(X) where [0, X] is the interval of integration (fractional 
differentiation) with 0 < x n < X. This algorithm was arrived at by approximating 
the integral (1.8c) with a trapezoidal product method, thereby restricting 0 < a < 2. 
Similar algorithms applicable to larger ranges of a can be constructed by using the 
general procedure derived in Ref. [20], if they become needed. 

The Griinwald-Letnikov algorithm is often used to numerically approximate the 
Riemann-Liouville fractional derivative (cf., e.g., with Oldham and Spanier [82, §8.2] 
and Podlubny [86, Chp. 7]) and it was the first algorithm to appear for approximating 
fractional derivatives (and integrals). 

The extent of rememberance of past states exhibited by the hereditary nature 
of a fractional derivative is manifest, for example, in its weights of quadrature, as 
illustrated in Fig. 1.1. This operator exhibits a fading memory: 0.001 < |a 8j8 | < 0.01 
for the six cases plotted in this figure. If Dy(X) were to be approximated by a 
backward difference with h — A/8, then the effective weights of quadrature would 
be a 0>8 = 1 and ai i8 = — 1 with all remaining weights being zero, as represented 
by the line segments in this figure. Similarly, if D 2 y(A) were to be approximated 
by a like backward-difference scheme, then o 0)8 = 1, oi, 8 = —2 and 02 , 8 = 1 with all 
remaining weights being zero. It is evident from the data presented in Fig. 1.1 that the 
weights of quadrature a„ ]8 for approximating D°y(X) are compatible with those for 
the first- and second-order backward differences, and that fractional quadratures have 
additional contributions that monotonically diminish with increasing nodal number 
from node n = 2 fading all the way back to the origin at node n = N. This suggests 
that a truncation scheme may be able to be used to enhance algorithmic efficiency 
for some classes of functions, but not all. 


§ Apparently this algorithm first appeared in the PhD thesis of Chern [15], unbeknownst to us 
(KD) at the time of writing Ref. [20]. Chern used this algorithm to differentiate a Kelvin- Voigt, 
fractional-order, viscoelastic, material model in a finite element code. He did not address stability 
or uniqueness of solution issues; he did not compute error estimates; and he did not utilize an 
extrapolation scheme to enhance solution accuracy. 
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Algorithm 1.1 Computa tion of a Caputo fractional derivative (0<a<2,a/l). 
For interval [0, X] with grid {x n = nh : n — 0, 1, 2, , N} where h = X/N, compute 

T~\a fu\ 1 ( V''L Q J ( N-n) k h k (fc)\ 

-D* y N{h) — h a T(2-a) £r»=0 a n,iV [YN-ti 2^fc= 0 fe! ^0+ J’ 

D?y(X) = D?y N (h) + 0(h?-°), 

using the quadrature weights (derived from a trapezoidal product rule) 

{ 1, if n = 0, 

2 1-Q - 2, if n = 1, 

(n + I) 1 ' 0 - 2n 1_a + (n - l) 1_a , if 2 < n < TV - 1, 

(1 - a)N~ a - N 1 -** + {N- l) 1 - 0 , if n = N. 

Refine, if desired, using Richardson extrapolation 

£>*y“ = (D?y u v ll - 2 r " _1 D? yr 1 ) / (1 - 2 r “-), 

D*y(X) = D“y“ + 0(h ru ), 

such that if 0 < a < 1 then r u _i is assigned as 

r 0 = 2- a, 
r\ = 2, r 2 = 3 - a, r 3 = 4 - a, 
r 4 = 4, r 5 = 5 — a, r 6 = 6 - o, 

r 7 — 6, 


Diethelm’s Quadrature Weights for Fractional Differentiation 

x /X, 0 < x < X 



Figure 1.1: Weights of quadrature a„ iN for approximating Caputo’s fractional deriva- 
tive (1.8) over interval [0,X] using Diethelm’s [20] Alg. 1.1, plotted here for various 
values of a with N = 8. 
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Richardson extrapolation 

Richardson extrapolation is a technique that can often be used to increase the accu- 
racy of results [24], As we utilize it, this technique follows a triangular scheme — a 
Romberg tableau — that has the form 


d? y 0 ° 
DS y? 
D?y° 2 
D?y° 3 


D *y\ 
D?y l 
D?y l 


D?yl 

D *yi 


D?yl 


( 1 . 10 ) 


Constructing the first column of the tableau constitutes the bulk of the computational 
effort. Here D* y° := D/y^(h) denotes the value of D°y evaluated numerically at X 
over [0, AC] using an initial stepsize of h (= x /n ) , D* y° := D+yu( h / 2 ) is computed 
using the refined stepsize of ~h (= x / 2 n), D“ y° := £>“ y,/v(V4) is computed using the 
more refined stepsize of \h (= x /an), while := yN( h /&) is computed using 

the further refined stepsize of |/i (= x /&n), etc. The remaining columns are quickly 
evaluated using the recursive formula listed in Alg. 1.1. The advantage of constructing 
this tableau is that the elements D° y“(X) in the u th column converge for fixed u and 
increasing v towards the true value of the Caputo derivative as 0(h ru ). Hence, the 
further one moves to the right in the tableau the faster the column converges, and 
this level of convergence requires less computational effort to achieve than a direct 
computation of D°y N (X; h) when computed to a similar accuracy of 0(h r °) ~ 0(h ru ). 

Step-size choice 

The error analysis mentioned above is only a truncation error analysis. It assumes 
that the calculations are done in exact arithmetic, and it does not take into account 
effects like roundoff. When one needs to look at these effects too, it is possible to 
ask for a step size h = X/N whose combined effect arising from both error sources is 
minimized. As we have seen above, it is likely that the truncation error decreases with 
the step size h, whereas roundoff tends to have the opposite behavior, so we should be 
looking for a sort of compromise. The considerations in this context are very similar 
to those for integer-order derivatives [89, §5.7]. Roughly speaking, it turns out that 
the roundoff error behaves as h~ a e y f(r]), where e y is the relative accuracy with which 
one can compute y, and where 77 designates some number within the interval [0, X], 
Moreover, the truncation error is close to c a h 2 ~ a f'( £), where c Q is some constant 
independent of /, and £ is some other number also contained in the interval [0,X]. 
Consequently, an optimum step size would be of order h ~ when 

minimizing with respect to both trucation and roundoff errors. 

Unless specific information indicating the contrary is available, one may assume 
that / and f" are not too irregularly behaved. Under these conditions /(r/) ~ f{X) 
and f'{r]) « f"(X), and one can then follow the suggestion of Press et al. [89, p. 187] 
by setting f(X)/f"(X) « X (except near X = 0 where some other estimate for this 
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quantity should be used). This scheme provides some advice on the choice of step size 
if roundoff effects are considered problematic in some specific application at hand. 

1.4.2 Riemann-Liouville Fractional Integrals 

In the course of our work we shall not only have to approximate fractional derivatives, 
but also fractional integrals. As indicated above, the natural concept for the fractional 
integral to be used in connection with Caputo derivatives is the Riemann-Liouville 
integral described in (1.2). We therefore present a numerical scheme for the solution 
of this problem, too. The underlying idea of the algorithm, stated in a formal way 
in Alg. 1.2 below, is completely identical to the idea presented above for the Caputo 
derivative; that is, we use a product integration technique based on the trapezoidal 
quadrature rule. Said differently, we replace the given function f on the right-hand 
side of (1.2) by a piecewise linear interpolant, and then we calculate the resulting 
integral exactly. As a matter of fact, this algorithm will also be part of the scheme 
introduced in Alg. 1.3 in the pages that follow for the numerical solution of certain 
Caputo- type differential equations. 

It is easily seen that the error of this algorithm is of the order 0(h 2 ) where, as 
above, h denotes the step size. Once again, we can improve the accuracy by adding a 
Richardson extrapolation procedure to the plain algorithm. The required exponents 
are known (cf. [52, §4]) and the resulting scheme is detailed in Alg. 1.2. Both the 
fundamental algorithm itself, and the Richardson extrapolation procedure, may be 
used for any positive value of a ; there is no need to impose an upper bound on the legal 
range for a. This is due to the fact that the Abel (power law) kernel in the definition 
(1.2) of the Riemann-Liouville integral is regular, or at worst, weakly singular, and 
hence, integrable (at least in the improper sense) for any a > 0. In contrast, the 
corresponding kernel in the definition (1.8c) of the Caputo derivative is not integrable. 
This kernel requires special regularization methods that are compatible with our 
approximation method, and as such, our scheme for approximating Caputo derivatives 
is only valid for 0 < a < 2, whereas, our corresponding scheme for approximating 
Riemann-Liouville integrals is valid for all a > 0. 

Notice the formal correspondence between Alg. 1.2 (for fractional integration of 
order a) and Alg. 1.1 (for fractional differentiation of order a). Except for the 
initial conditions that have to be taken into account additionally, the latter is simply 
obtained from the former by replacing the parameter a by — a.^ This relates to 
the intuitive (but not mathematically strictly correct) interpretation of fractional 
differentiation and integration being inverse operations. Also notice that the index 
ordering is inverted between these two algorithms, which is in keeping with accepted 
indexing conventions. Algorithm 1.1 indexes from xo = X to = 0, while Alg. 1.2 
indexes from xo = 0 to = X. 

A visualization of quadrature weight versus nodal index for several values of a 
pertaining to Alg. 1.2 is presented in Fig. 1.2. What is striking about this figure is the 

^Similarly, the Griinwald-Letnikov algorithm for approximating Riemann-Liouville fractional 
derivatives of order a also applies for approximating Riemann-Liouville fractional integrals by re- 
placing their algorithmic parameter a with —a [82, §8.2]. 
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Algorithm 1.2 Computation of a Riemann-Liouville fractional integral (a > 0). 

For interval [0, A] with grid {x n = nh: n = 0, 1, 2, , N} where h — X/N, compute 
J a yN(h) = r(2 +Q ) En=0 C n,N Yn, 

J a y(X) = J a y N (h) + 0(h 2 ), 

using the quadrature weights (derived from a trapezoidal product rule) 
f (1 + a)N Q - N l+a + (N — l) 1+a , if n = 0, 

c n>N = < (N-n+l) 1+a ~2(N-n) 1+a + (N-n-l) 1+a , if 0 < n < N, 

^ 1, if n = N. 

Refine, if desired, using Richardson extrapolation 

J a r v = - 2 r “~ l j° yr 1 ) / (1 - 2 r “- 1 ). 


J a y(X) = J a y u u + 0(h r »), 
such that if 0 < ct < 1 then r u _i is assigned as 


-5 

o 

II 

to 

r i 

= 2 + a 

r 2 = 3, 

r 3 

= 3 + Oi 

f 4 = 4, 

r 5 

= 4 + a 

r 6 = 5, 




Note: Whenever a > 1, the same values appear in the sequence r 0 , r\, r 2 , . . ., but 
they now have to be ordered in a different way to keep the sequence monotonic. (For 
example, if 1 < a < 2 then we have r 0 = 2, 7 ^ = 3, r 2 = 2 + a, r 3 = 4, r 4 = 3 + cm , ■ • -)• 


obvious difference between domains 0 < a < 1 and 1 < a < 2. Whenever a = 1, the 
algorithm reduces to classic trapezoidal integration. Whenever 0 < a < 1, the earlier 
states will contribute less to the overall solution than will the more recent states, but 
they do not entirely fade out. Fractional integration exhibits long-term memory loss 
when 0 < a < 1 but, unlike fractional differentiation, fractional integration does not 
experience a total loss (or fading away) of past memories. Also, the smaller the value 
of a (i.e., the closer it is to zero) the greater the degree of long-term memory loss 
will be. In contrast, whenever 1 < a < 2, the earlier states will contribute more to 
the overall solution than will the more recent states. Fractional integration therefore 
exhibits short-term memory loss when 1 < a < 2. This is like an elderly person who 
remembers in vivid detail what happened years ago, but who cannot recall what took 
place yesterday. Furthermore, the greater the value of a (i.e., the closer it is to two) 
the more pronounced the short-term memory loss will be. 

The line segments displayed in Fig. 1.2 represent averaged and normalized weights 
of quadrature over each subinterval. The actual nodal weights, c H} n, are often ob- 
served to be non-monotonic at either of the two nodal endpoints. In this integration 
scheme there are N - 1-1 nodal weights that apply to N subintervals, but there should 
be exactly one weight per subinterval. So how the algorithm works (internally, and 
roughly speaking) is to average these weights in a trapezoidal fashion, as outlined in 
Table 1.1. In other words, the inner weights are divided into two equal halves with 
each half going to one of the two adjoining subintervals. In addition to averaging, the 
displayed line segments have been normalized to the interval [0, 1]. Normalizing allows 
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Diethelm’s Quadrature Weights for Fractional Integration 

x/X, 0 < x < X 

0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 



Figure 1.2: Effective weights of quadrature (c n>N ) for approximating the Riemann- 
Liouville fractional integral (1.2) over interval [0, X] using Alg. 1.2, plotted here for 
various values of a with N = 8. 


Subinterval 

Averaged Quadrature Weight (c„,jv) 

[0,X/N\ 

(co.jv) = C 0 ,jV + I Ci.jV 

[nX/N, (n + 1 )X/N\ 

(Cn,iv) 2 (pn,N 1) 2, ... , N 2) 

[(N-l)X/N,X] 

(CN-1,n) = \ Cjv— 1, N + CAT, AT 


Table 1.1: Averaging procedure used to compute effective weights of quadrature for 
approximating Riemann-Liouville integration as they relate to Alg. 1.2. 


one to discern the influence of a on quadrature in a meaningful way. The outcome is an 
averaged and normalized quadrature weighting that is monotonic in the nodal index 
number, as demonstrated by the line segments in Fig. 1.2, where there is a monotonic 
increase (decrease) in the effective weight of quadrature, (c n> jv) / max m (c TO! Ar), with 
increasing nodal index number for 0 < a < 1 (1 < a < 2). 

1.4.3 Caputo-Type FDE’s 

A numerical algorithm that solves Caputo-type FDE’s has been derived by Diethelm 
et al. [23] and is listed in Alg. 1.3. A thorough analysis of its algorithmic error is given 
in [22], This algorithm is of the PECE (Predict-Evaluate-Correct-Evaluate) type. 
Other numerical algorithms exist that solve FDE’s (e.g., Gorenflo [44] and Podlubny 
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Algorithm 1.3 Computation of a Caputo FDE (0 < a < 2, a ^ 1). 

For interval [0, X] with grid {x n = nh : n = 0, 1, 2, . . . , N : h = X/N}, predict with 

y£0) = Eil & yj? + ftfcy ESf f(*„, y»), 

using the quadrature weights (derived from a rectangular product rule) 
b n , N = {N-n) a -{N-n- 1) Q , 
and evaluate f(X,y^), then correct with 


(jZnJ C ",N f (*», y») + C N ,AT f(X, y£)) , 


yjv(/») = Ei=o f yj+ + f(&y (E^o 1 c n ,jv f (*„, y») + f(X, y £) ) , 
y(X) = yiv(h) + 0(h min ( 1+a ’ 2 )) , 

using the quadrature weights (derived from a trapezoidal product rule) 

' (l + a)N a - N 1+a + (N -l) 1+a , if n = 0, 

c„ >JV = (A — n+l) 1+Q — 2(N — n) 1+a + (N —n — l) 1+a , if 0 < n < N, 
\l, if n = N, 

and finally re-evaluate f(X,y^) saving it as f(£jv,y.v)- 

Refine, if desired, using Richardson extrapolation 

y“ = (y^i 1 - 2 r “-‘ y r 1 ) / (i - 2,w )> 

y(A) = y“ + 0(/i r “), 

such that whenever 0 < a < 1 the exponent r u _i is assigned as 

r 0 = 1 + a, 
r x — 2, r 2 = 2 + a, r 3 = 3 + a, 

r 4 = 4, r 5 = 4 + a, r 6 = 5 + a, 

r 7 = 6, 

or whenever 1 < a < 2 it is assigned as 
r 0 = 2, ri = 1 + a, r 2 = 2 + a, 

r 3 = 4, r 4 = 3 + a, r 5 = 4 + a, 

r 6 = 6, 


[86, Chp. 8]), but they focus on solving Riemann-Liouville FDE’s and usually restrict 
the class of FDE’s to be linear with homogeneous initial conditions. Algorithm 1.3 
solves non-linear Caputo FDE’s with inhomogeneous initial conditions, if required. 

The restriction that a ^ 1 in this algorithm is purely for formal reasons. If 
a — 1, then we can still implement the algorithm exactly in the indicated way. It 
must be noted, however, that it then is the limit case of an algorithm for fractional 
differential equations, and these equations involve non-local differential operators. 
Thus, the resulting scheme is non-local, too. In contrast, a method constructed for a 
first-order equation will, in practice, always make explicit use of the local structure of 
such an equation to save memory and computing time. Therefore, the case a — 1 of 
our algorithm will never be a competitive alternative to the usual methods for first- 
order equations. In particular, our algorithm is distinct from the algorithm known as 
the second-order Adams-Bashforth-Moulton technique for first-order problems when 
a = 1 . 
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Diethelm’s Quadrature Weights for Fractional Integration 

x/X, 0<x<X 

0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 



Figure 1.3: Normalized weights of quadrature b Ut jv for the predictor that approximates 
Caputo FDE’s (1.9) over interval [0,X] when using Diethelm et al.’s [23] Alg. 1.3, 
plotted here for various values of a with N = 8. 


Illustrations of quadrature weight versus nodal index for several values of a, as 
they pertain to the PECE method of Alg. 1.3, are presented in Figs. 1.2 & 1.3 with the 
former figure pertaining to the corrector and the latter pertaining to the predictor. 
FDE’s, like fractional integrals, exhibit long-term memory loss when 0 < a < 1, no 
memory loss when a = 1, and short-term memory loss when 1 < a < 2. 

Unlike the (c n>N ) in Fig. 1.2, where the N + 1 quadrature weights are averaged 
at the beginning and end of each subinterval in order to get N effective weights 
for N subintervals, the b UtN in Fig. 1.3 are fixed to the beginning of each of the N 
subintervals, and as such, do not require any ‘effective averaging’ to take place. This 
is a consequence of the b n<N quadrature weights belonging to an explicit integrator, 
while the c n< jv weights belong to an implicit integrator. Contrasting Figs. 1.2 & 1.3, 
there is little difference between the b n> n and (c Ui n} curves, indicating that there is a 
much stronger influence of a on the weights of quadrature than there is on the order 
of accuracy (e.g., 0(/i mm ( 2,1+a ))) that a particular integration scheme provides. 

Differential equations of fractional order have found recent applications in a variety 
of fields in science and engineering (e.g., see references in [57, 86]): chemical kinetics 
theory, electromagnetic theory, transport (diffusion) theory, fractal theory, control 
theory, electronic circuit theory, porous media, etc. One of the first applications of 
the fractional calculus was viscoelasticity, which is the primary focus of this work. 
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Efficient approximations 

Ford and Simpson [35] have extended Alg. 1.3, which possesses 0(n) operation counts 
at each stage and 0(n 2 ) overall, to a more efficient scheme with 0(n log n) counts 
overall, while retaining the accuracy of the method. 


1.5 Mittag-Leffler Function 

The (generalized) Mittag-Leffler function E a ^(z) is an entire function (in z G C) of 
order 1 / a that is defined by the power series [31, §18.1] 

00 k 

S..o(0 = S r w Z +ak y < 1U ) 

with E a (x) = E a< i(x) being the original function studied by Mittag-Leffler [79]. This 
function plays the same role in differential equations of fractional order as the expo- 
nential function e z plays in ordinary differential equations; in fact, Ei t i(z) — e z . 

A special form of the Mittag-Leffler function, 

:=£ q> i (-{(*- 0/r)“), 0 < a < 1, 0 < r, 0 <t'<t, 

and its derivative, 

m - f) := 8G( * ~ ^ B „, 0 (- ((1 - 

appear in FOV, which is the subject of much of this report. 

We now present some important properties of the Mittag-Leffler function, and a 
numerical algorithm for its rigorous solution, both of which are useful when consid- 
ering differential equations of fractional order. 

1.5.1 Analytical Properties 

In spite of the fact that in applications to differential equations of fractional order 
where the Mittag-Leffler function is typically restricted to the real line, we still need 
to give some of its properties in the complex plane. The main reason for this is that 
the numerical algorithm we present in the next sub-section consists of two parts: the 
first part gives a numerical value for the Mittag-Leffler function with a < 1, while the 
second one uses, for a > 1, some special formulae that reduce this case to the previous 
one. These special formulae are defined over the complex plane and are given by 

1 m 

E ° Az) = 2^TT £ "> = 0,1,2,..., (1.12a) 

h=—m 

and 

1 m — 1 

1 £ E a /m ,p(z 1/m e 2< " h/m ) , m — 1,2, , (1.12b) 
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where i (:= \/— l) is the imaginary unit number. Obviously, if a > 1, and even if z 
is a real number, we still need to evaluate the numerical values of the Mittag-Leffler 
function with, in general, a complex argument (and a < 1). 

First, we present some important integral representations of the Mittag-Leffler 
function. Let us denote by 7 (e; p) (e > 0, 0 < p < 7r) a contour in the complex 
A-plane with non-decreasing arg A consisting of the following parts: 

1) the ray arg A = — p, |A| > e; 

2) the arc — p < arg A < p from the circumference |A| = e; 

3) the ray arg A = p , |A| > e. 

In the case where 0 < < 7r, the complex A-plane is divided into two unbounded parts 

by the contour 7 (e; p): domain p) is to the left of the contour, while domain 

GW(e; p) is to its right. If p = n, the contour y(e; p) consists of the circumference 
|A| = e and the cut — oo < A < — e. In this case, domain p) becomes the circle 

|A| < e, while domain G^ + ^(e; p) becomes the region {A : I arg A| < 7 r, |A| > e}. 

Let 0 < a < 2, let /3 be an arbitrary (real or complex) number, and let a non- 
negative number p be chosen such that 

an 

— < p < mm{7r, an). (113) 

z 

Then we have the following integral representations for the Mittag-Leffler function: 


B ° dz) ~ 2 iria L 






A — z 


dX, z 6 \e\p), (1.14a) 


and 


„ , , _ z<-»/.e' v ” 1 

Ea,p{z) — b 


a 


— f 

2 nia J yU , 




A — z 


dX, z € G^ + \e\ p) . (1.14b) 


If j5 is a real number, as assigned in (1.11), then the formulae of (1.14) can be 
rewritten in forms that are more suitable for numerical evaluation (see Gorenflo et al. 
[45]). In particular, if 0 < a < 1, /3 € R, j arg z| > an, z / 0, then 


EaA Z ) = 

J 

p OO pOCK 

/ K(a,p,x,z)dx+ P(a,p,e,p,z)dp, e > OJel, 

J —an 

(1.15a) 


poo 

E at/3 (z)= K(a,fi,x,z)dx, if /3 < 1 + a, 

Jo 

(1.15b) 

Ea,0 (^) — 

sin(aTr) f°° e~x l,a J 1 . f 0 , , 

an J 0 X 2 ~2 X zcos{an) + z 2 X z’ 1 ^ + 

(1.15c) 


where 


K(a,(3,x,z ) 


X (1 0) G e x 1/a ^sin(7r(l —/?))— z sin(7r(l — (3 + a)) 
an x 2 ~ 2x z cos(an) + z 2 


P(a,/3,e,p,z ) 


e 1+(1 0)/a e e 1/Q cos(»ya)( cos ^ + isin(o;)) 
2o:7r ee iip — z 
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oj = e Vo sm(v/a) + ip(l + ( 1_ ^)/a). 

The representations in (1.15), and similar formulae for the cases | arg z\ = air and 
| arg z\ < an presented in Gorenflo et al. [45], are an essential part of the numerical 
algorithm listed in the next sub-section. 

Using the integral representations in (114), it is not too difficult to get asymptotic 
expansions for the Mittag-Leffler function in the complex plane. Let a < 2, f3 be an 
arbitrary number, and <p be choosen to satisfy the condition (1.13). Then we have, 
for any p € N and \z\ — > oo: 

1) Whenever | argz| < <p, 




,<1 -V/apZ 1 '* 


a 


v z 




(1.16a) 


2) Whenever <p < | argz| < 7r, 


B ^w = -Er(y^) + 0 W I_P )- ( 116b > 

These formulas are also used in our numerical algorithm. 

Thorough discussions of properties of the Mittag-Leffler function can be found, 
for example, in Refs. [31, 76, 86]. 


1.5.2 Numerical Algorithms 

Robust 

The numerical scheme listed in Alg. 1.4 for computing the general Mittag-Leffler func- 
tion E a ^(z) is taken from an obscure paper written by Gorenflo et al. [45]. Their 
algorithm uses the defining series (Eqn. 1.11) for arguments of small magnitude, its 
asymptotic representation (Eqn. 1.16) for arguments of large magnitude, and special 
integral representations (the formulae in (1.15) for the case where | arg z\ > an, and 
similar representations for the cases | arg z\ = an and | arg z | < an) for intermediate 
values of the argument that include a monotonic part f K(a, 13, x, z) dx and an os- 
cillatory part f P(a, j3, e, <f>, z) dcj), which can themselves be evaluated using standard 
techniques (cf. App. B). 

Efficient 

Algorithm 1.4 can produce a numerical result to any desired level of accuracy, but 
these computations are expensive and therefore their use in a finite element setting, for 
example, is prohibative. To meet this need, we have constructed a table of Pade ap- 
proximates for E a (—x a ) in App. C for x > 0 and a 6 {0.01, 0.02, 0.03, . . . , 0.98, 0.99}. 
As we shall see in the next chapter, this form of the Mittag-Leffler function arises in 
many fractional-order, viscoelastic, material models, including those of interest to us. 

Another algorithm for solving the Mittag-Leffler function E a (x) (0.02 < a < 0.98 
with a reported relative error that is less than 1.6 x 10 -5 ) has been published by 
Welch et al. [109]. 
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Algori thm 1.4 Computation of the Mittag-Leffler function. 

GIVEN a > 0, p eR, zeC THEN 
IF 2 = 0 THEN 

E at p(z) = ppy 

ELSIF \z\ < 1 THEN 

fc 0 = max{ f(i-0)/o"| , fln[e m (l- \z\)\/ ln(|z|)] } 

E a ,p(z) = Ylk=0 r(P+ak) 0( £ m) 

ELSIF \z\ > [10 + 5ctJ THEN 
k 0 = [— ln(e m )/ ln(|z|)J 
IF | argz| < aiT /b + V 2 min{7r, an} THEN 

- T.IU fwUm + 

ELSE 

E a fi{z) — — Ylk= 1 T(J}-ak) 

ELSIF a < 1 THEN 

_ ( max{l, 2\z\, (-ln(* e »»/6)) a }, P > 0 

X ° _ 1 max{(|/3| + l) a , 2\z\, (—2 ln(^-/ [6( |^| +2)(2 |>9| ) i>®!] ))“ }, P < 0 
K(a,P, X ,z) = 

P(a,P,eA,z) = exp( e V« C os(Va))^g^ 

U = 4 >( 1 + ( 1_ ^)/a) + el/ “ Sin(^/ a ) 

IF | argz| > an THEN 
IF p < 1 THEN 

E a #(z) = f 0 Xo K(a, p, x, z) dx + 0(e m ) 

ELSE 

E a ,fs{z) = f X0 K(a, p, x, z) dx + P(a, P, 1, <t>, z) d(f> + 0{e m ) 

ELSIF | argz| < an THEN 
IF p < 1 THEN 

Ea,p(z) = f xo K(a,p,x,z)dx+ ^z {1 ~ e) ^e zl/a + 0{s m ) 

ELSE 

£«*(*) = C /2 /?, X, *) ^X + Hr & |2| /2^> 2 ) d( t> + i* <1H,)/ac * 1 

+ G(e m ) 

ELSE 

£«*(*) = iS +1)/ , *(«, P, x , z) d X + /_°1 P(a, P, (W+D/ 2 , <t > , z) d* + 0(e m ) 
ELSIF 1 < a < 2 THEN 

-£'a,/3(- 2 0 = Ea^fi^Z 1 ! 2 ) + E a /2fi{ — Z^ 2 ) 

ELSE 

= [ a h\ + 1 

SaJiM = E Etc' «./*.* (* V *» “P( 2 ’“A.» 

END 

Parameter e TO denotes machine epsilon (precision). 
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Chapter 2 
ID FOV 


In the 1940’s, Scott Blair [8] and Gerasimov [42] independently proposed a material 
model bounded between a Hookean solid (a = 0) and a Newtonian fluid (a = 1). 
Their relationship — a fractional Newton model — can be written as <r(t) = [xT a D°s(t), 
where a and e denote stress and strain, respectively, which are considered here to be 
causal functions of time t. The coefficient jut (> 0) represents a single material 
constant (a generalized viscosity: /j, has units of stress, while r has units of time), 
and exponent ct (0 < g: < 1) can be considered as a second material constant. 
Experimental results motivated Scott Blair’s model development. Mathematics, on 
the other hand, motivated Gerasimov who was the first to consider an Abel kernel 
for the relaxation modulus in Boltzmann’s general theory of viscoelasticity. 

Bagley and Torvik [4] demonstrated that the molecular theory of Rouse (for dilute 
solutions of non-crosslinked polymer molecules residing in Newtonian solvents) has 
a polymer contribution to stress that corresponds to a fractional Newton element 
whose order of evolution is a half (i.e., a = ^/i)- They also state (without proof) that 
the molecular theory of Zimm (for dilute solutions of crosslinked polymer molecules 
residing in Newtonian solvents) has a polymer contribution to stress that corresponds 
to a fractional Newton element whose order of evolution is two thirds (i.e., ot = 2 /s)- 

Gemant [38] was the first to propose a fractional viscoelastic model. He extended 
the notion of a Maxwell fluid by replacing its first-order derivative on stress with the 
semi-derivative, and in doing so, he proposed that [l + \J 77/ \xDJ ]c(t) = r]D£(t), 
where p, (> 0) and 77 (> 0) are material constants. The fractional Maxwell fluid, 
which is a spring in series with a fractional Newton element, actually has the form 

[1 + T a D°]a(t) = T ] T a - 1 D«e{t), <7 0+ = ^£ 0 +, (2-1) 

where 77 (> 0) is the viscosity, r (> 0) is the characteristic relaxation time, and 
exponent a (0 < a < 1) is the fractal order of evolution, which is taken to be the 
same for both stress and strain, while <tq+ and £0+ are 4he initial states of stress 
and strain at time t = 0 + , thereby allowing for an inhomogeneous initial state of 
fin ite stress — a characteristic that Gemant’s model does not possess. The fractional 
Maxwell fluid was first discussed in the manuscript of Caputo and Mainardi [13] as 
a special case to their material model (Eqn. 2.2 below). We refer to (2.1) as the 
standard FOV fluid in ID. 
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Caputo [12] introduced a fractional Voigt solid a(t) = /j[l + p a D^]£(t) to model 
the nearly rate-insensitive dynamic response of Earth’s crust over large ranges in 
frequency when excited by earthquakes. Here p (> 0), p (> 0) and a (0 < a < 1) 
are the material constants. As a mechanical model, this is a spring in parallel with 
a fractional Newton element. A more appropriate representation of solid behavior is 
the fractional Kelvin model, which is a spring in parallel with a fractional Maxwell 
element. This material model was introduced by Caputo and Mainardi [13] and has 
the form 

[l + T a D°]a(t) = p[l +p a D“]s(t), a 0 + = p £ 0 +, (2.2) 

where p (> 0) is the rubbery modulus, /x(p/r) Q (> p) is the glassy modulus, r (> 0) 
is the characteristic relaxation time, p (> r) is the characteristic retardation time, 
and exponent a (0 < a < 1) is the fractal order of evolution. This model, unlike 
Caputo’s original model, allows for an inhomogeneous initial state of finite stress. 
Bagley and Torvik [5] have shown that the fractal orders of evolution in stress and 
strain must be the same, as written in (2.2), and as originally proposed by Caputo and 
Mainardi, in order for this constitutive realationship to be compatible with the second 
law of thermodynamics; specifically, in order to guarantee a non-negative dissipation 
whenever a cyclic loading history is imposed on the material. We refer to (2.2) as the 
standard FOV solid in ID, in the spirit of Zener [111, pg. 43] referring to Kelvin’s 
model [1 + rD\a{t) = p[ 1 + pD]£(t) as the “standard linear solid”. 

The initial conditions present in (2.1 & 2.2) come about by taking the Laplace 
transform* of these constitutive formulae. What one learns from these transformations 
is that if the material model is to be physically admissible, in the sense that it 
propagates a wave front at finite speed, then that part of the transformation which 
pertains to the initial state must be independent of the Laplace transform variable s in 
the frequency domain, or it must have like dependencies on both sides of the equation 

*The Laplace transform f(s) of function f(t) is given by the mapping procedure f(t) 4 /(.s) = 
f 0 exp(— st) /(t) dr, where 4 denotes the juxtaposition of function /(f) with its Laplace transform 
/(s). In fractional-order viscoelasticity, the Laplace transform pairs 


[aj 




k = 0 


t 71 " 1 ^ J_ 
T(n) ’ s n 


and t 13 1 E a! p(±at a ) 


S a-P 
s a =F a 


have particular significance, where a, a,n, t € R+, ft € R and E a ,p(t) = J2T=o /E(/? 4 ak) is the 
general Mittag-Leffler function, which plays a role in FDE’s like that which the exponential function 
plays in ODE’s; in fact, Ei,i(t) = e‘. 

The above formulae are analytic continuations of the well known Laplace transform pairs 




-7, (fc) 


*771 — 1 


0 + ’ 


k=0 


[m — 1)! 


4- — and 
s m 


e ±at 4 


1 

S -p Qi 


where m G N and a, t e 1+ . In contrast to the Laplace transform of Caputo’s deriviative, which 
contains a sum of integer-order derivatives of the initial state, the Laplace transform of the Riemann- 
Liouville derivative contains a sum of fractional-order derivatives of the initial state, making the 
initialization of Riemann-Liouville based differential equations a difficult task, but not an impossible 
one [73]. 
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that then cancel out in the initial state. Having derivatives of equal order on both sides 
of the equation, as in the standard FOV fluid and solid material models, is one way 
to ensure that this physical constraint is adhered to. In the aerodynamics literature, 
this process of addressing the initial state for consistency of initial condition in the 
frequency domain is known as the method of shocks, which was introduced into the 
FOV literature by Bagley and Calico [3]. Another very important reason to restrict 
the class of admissible material models to only include those that propagate waves 
at finite speeds has to do with stability. Material models that predict infinite wave 
speeds will become mathematically unstable at some critical finite velocity [65]. 

One objective of this paper is to derive 3D versions of the standard, FOV, fluid 
and solid, material models without imposing any constraints as to the magnitude of 
deformation. To the best of our knowledge, Drozdov [27] is the only person to have 
extended linear, fractional-order, viscoelastic formulations into 3D formulae applicable 
to non-linear mechanics where finite strains are present. Specifically, he extended 
the following ID models: [l + {r]/p) a D^]a(t) — rjDs(t) y which is a generalization 
of Gemant’s [38] fractional Maxwell model, and a(t) = /u[l + p a D^]e(t), which is 
Caputo’s [12] fractional Voigt model. In Chps. 7-??, we introduce 3D versions for 
the standard FOV fluid and the standard FOV solid, which are presented here in ID 
in Eqns. (2.1 & 2.2). 


2.1 Material Functions 


The parameterization procedures that follow assume infinitesimal strains in homoge- 
neous ID deformations. 

Boltzmann’s [9] linear theory of viscoelasticity, which includes the standard FOV 
models of (2.1 & 2.2), can be expressed as an integral equation with a hereditary 
kernel that convolves with a change in the independent state variable according to 
the convolution rules of either Stieltjes or Riemann. Whenever stress responds to 
strain, this theory can be expressed in terms of a (relaxation) modulus G(t) where 
[16, pp. 3-9] 


a (t) = f G(t - t') ds(t') = e 0 +G(t) + [ G{t - t ') De(t') dt', (2.3a) 

Jo J o+ 

or conversely, whenever strain responds to stress, Boltzmann’s theory can be re- 
expressed in terms of a (creep) compliance J(t ) where 

e(t) = r J(t - t') dcr(t') = (Jo + J(t ) + [ J(t- t') Da(t') dt'. (2.3b) 

Jo J o+ 

These two convolution integrals can be solved analytically using Laplace transform 
techniques, provided that the loading histories are simple enough. 

The standard FOV fluid (2.1) has a modulus G(t) and a compliance J(t) of [13] 


G(f) =?£„(-(</,)”) 


(*A)° 

r(l + a) 


Jit) — — 1 + 




(2.4) 
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whereas, the standard FOV solid (2.2) has the material functions [13] 




i i+ r 


p° — T a 


where E a (x) = E a<1 (x) is the Mittag-Leffler function (see §1.5). 


(2.5) 


2.1.1 Static Experiments 

Stress relaxation experiments are often executed for the purpose of materials charac- 
terization, where 

a(t) j E a [— (Vr) a ) for standard FOV fluids, 

0o+ [^ + pa ~J a E a (—( t / r ) a ) for standard FOV solids. 

Figure 2.1 presents a normalized plot of stress relaxation curves for the standard 
FOV fluid, with a = 1 designating the response of a classic Maxwell fluid. The 
stress relaxes to zero monotonically in a FOV fluid for all a € (0,1]. Figure 2.2 
presents a normalized plot of stress relaxation curves for the standard FOV solid, 
with a = 1 designating the response of a classic Kelvin solid. For all a e (0, 1], the 
stress monotonically relaxes to a unique non-zero value in a FOV solid as t — > oo, 
which distinguishes solid behavior from fluid behavior. Here, and in the following 
figures of this chapter, the relaxation r and retardation p times are assumed to scale 
as ( p/t ) q = 5 for purposes of illustration. These figures show that the fractal order of 
evolution controls the shape of the relaxation curve. 

Relaxation, as described above, exhibits an exponential decay as t — > oo when- 
ever a = 1, and an algebraic decay to infinity whenever 0 < a < 1. This corresponds 
to a regular rate process leading to strong mixing (exponential decay) versus an in- 
termittent rate process causing weak mixing (algebraic decay), as quantified by a 
probabilistic fluctuation of recurrent events (molecular collisions) governing the ve- 
locity relaxation process in polymer chain physics. Douglas [25] has shown, through 
probabilistic reasoning using Feller’s renewal theory, that the autocorrelation function 
describing relaxation phenomena is governed by a fractional-order differential equa- 
tion whose solution is given in terms of the Mittag-Leffler function, and whose order 
of evolution correlates with the degree of intermittency in the relaxation process. 

Douglas [25] also states that the stretched exponential, e.g., exp(— ( 4 / r ) Q ) , often 
used to empirically fit relaxation data in the literature, does not arise from probabilis- 
tic considerations in polymer chain physics. Popularity of the stretched exponential 
over the Mittag-Leffler function has two likely sources: many researchers are not fa- 
miliar with, or have even heard of, the Mittag-Leffler function, and if they are familiar 
with it, they do not likely know how to compute its value. With respect to the latter, 
see Alg. 1.4 and App. C. 
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Relaxation Modulus: Standard FOV Fluid 

[l +T “D>(t) = ^’0:6(1) 



Figure 2.1: Normalized diagram for stress relaxation of a fractional Maxwell fluid. 


Relaxation Modulus: Standard FOV Solid 

[1 + x“D“]a(t) = (x| 1 + p a D“]e(t) 



Figure 2.2: Normalized diagram for stress relaxation of a fractional Kelvin solid. 
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Creep Compliance: Standard FOV Fluid 

[1 + T a D“]<j(t) = Tyc a-1 D“e(t) 



Figure 2.3: Normalized diagram for creep of a fractional Maxwell fluid. 


Creep experiments are also performed for purposes of materials characterization, 
where 


g(0 

£o+ 


1 + ( t / T ) Q for standard FOV fluids, 

1 + p ^r- (l - E a (-0/ P ) Q )) for standard FOV solids. 


(2.7) 


Figure 2.3 shows a normalized plot of creep curves for the standard FOV fluid. A 
specimen will creep without bound in a FOV fluid for all a G (0, 1]. However, only 
in the case of a Maxwell fluid where a = 1 is the ‘effective’ viscosity (i.e., the slope) 
a constant. Conversely, FOV fluids will eventually (at infinite time) stop creeping 
altogether (see Eqn. 2.9). Figure 2.4 shows a normalized plot of creep curves for the 
standard FOV solid. Here creep stops at a unique threshold level in strain for all 
a G (0,1]. Steady-state creep behavior cannot be predicted by this class of material 
models. The fractal orders of evolution influence the shape of these curves, too, and 
in the case of a solid, they also influence the time required to attain saturation. 

It is difficult to parameterize an FOV solid with only relaxation data, or with only 
creep data. This is because it is difficult to acquire sufficient sensitivity in the data 
to the parameter p in the case of relaxation, or to the paramter r in the case of creep. 
But whenever relaxation and creep data are used together during estimation, ample 
sensitivty will exist for all material constants and good data fits can be expected. 

Although Figs. 2. 1-2. 4 are informative, they are not as practical as one would like 
in the sense that one cannot directly extract the order of evolution, a, from them via 
some graphical technique. However, if one were to measure stress rates in a relaxation 
experiment, or strain rates in a creep experiment, then the order of evolution could, 
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Creep Compliance: Standard FOV Solid 

[1 + tXWO = nil + p“D“]e(t) 



Figure 2.4: Normalized diagram for creep of a fractional Kelvin solid. 


in theory, be extracted as a slope in an appropriate log-log plot of the data. In a 
stress relaxation experiment 


Da(t) 
o o+ 


DE a ( ( */r ) a ) 

^ DE a (-(t/r) a ) 


for standard FOV fluids, 
for standard FOV solids, 


whereas, in a creep experiment 


De(t) 

£q+ 


TTifeTt'/.) 1 - 


for standard FOV fluids, 
for standard FOV solids. 


(2.8) 


(2.9) 


Figure 2.5 presents a graphical representation of dE a (—x a )/dx^ which appears in 
three of the four descriptions above, with the exception being creep rate in a standard 
FOV fluid, which has a power-law response. Cureously, what is observed in Fig. 2.5 is 
that dE a (-x a )/dx approximates power-law behavior whenever x < 0.1. The scheme 
depicted in this figure for extracting a is accurate (to two significant figures) over 
the range of 1 < a < x /i , but it looses accuracy as a approaches zero; for example, 
this graphical scheme yielded a value for a of 0.27 when it was actually Y 4 . Even so, 

^From Podlubny [86, pp. 21-22], one finds that 

dEg,l(-{*/y) a ) E a ,o{-(*/y) a ) 

dx x ’ 

for0<a<l,a:>0 and y > 0. 
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Derivative of the Mittag-Leffler Function 



Figure 2.5: Diagram of the derivative of the Mittag-Leffler function, dE a< i(—x a )/dx = 
E a fi{-x a )/x. 


this is likely to be accurate enough for material characterization purposes given the 
uncertainty of experiemental error. Because the time constant for creep is typically 
many orders of magnitude greater than the time constant for stress relaxation, it will 
be easier to satisfy a boundary of x < 0.1 in a creep experiment than it would be to 
satisfy it in a relaxation experiment; nevertheless, the experimental challenge remains 
great. 


2.1.2 Dynamic Experiments 

Also useful for the purpose of materials characterization are dynamic experiments 
where strain is controlled at a constant amplitude £ 0 and angular frequency u> ac- 
cording to e(t) = £ 0 exp(iivt), to which stress responds with a dynamic modulus of 
G*(oj) = G'(uj)-\-iG"(u ) such that a{t) = £qG*(uj ) exp {iust). The real G'{ u) and imag- 
inary G"(lo) parts of the dynamic modulus are called the storage and loss moduli, 
respectively, whose ratio, tan5(cu) = g "(")/g'(w) , is often reported in the literature. 
Figure 2.6 illustrates how these properties are extracted from experimental data ob- 
tained under steady, oscillatory, loading conditions, where the stress-strain curve is an 
ellipse with control e(t) = £ 0 sin cat and response a(t) = a 0 sin(cjf + tf). A material is 
non-linear if the hysteresis is something other than an ellipse under sinusoidal loading 
conditions. A thorough discussion of dynamic experiments, as they relate to linear 
viscoelastic materials characterization, can be found in the recent text by Lakes [62, 
Chp. 3]. 
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CT 



Figure 2.6: Extracting dynamic properties from linear viscoelastic materials. 


For the standard FOV fluid (2.1), the dynamic modulus is given by 

G ^ Tl + (iwr)“’ 


whose real and imaginary parts are 


V 


(wr) Q + cos ( an / 2 ) 


G ^ t ( wr)“ + (wr)-° + 2 cos( Q7r /2) 

r „ { , = V sin ( Q7r / 2 ) 

r (wr)“ + (o;r) _a + 2cos( a 7 r /2) 


and that ratio as 


tan5(w) = 


s in(aw/ 2 ) 


(wr)“ + COs( Q7r / 2 ) 

For the standard FOV solid (2.2), the dynamic modulus is given by 

-**/, _ .. (P\ a ( T / P ) a + (iur) a 


G 


r M = p(^) 


1 + (i u>t) 

whose real and imaginary parts are 

r' ( ( P Y ( UT ^ a + ^~ a + l 1 + ( T /^) Q ) C 0 S ( a V 2 ) "1 

\ UJ )-P\ T ) ( W r) a + (wT)-“ + 2cos(“V2) 

rHt \ = „ m a (i - (Vp) q ) sin(^y 2 ) 

^ ^ \t ) (cur) Q + (wr)" Q + 2 cos( a ’ 


< a V2) 


(2.10a) 


(2.10b) 


(2.10c) 


(2.11a) 


(2.11b) 
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Cole-Cole Plot: Standard FOV Fluid 

[1 + -c a D?]a(t) = r|x a ’ D“e(t) 



Figure 2.7: Normalized Cole-Cole diagram for a fractional Maxwell fluid. 


and that ratio as 


tan5(w) 


(l - i T / P ) a ) sin(°y 2 ) 

(WT ) Q + (up)~ a + (l + ( T /p ) a ) COs( a7r /2) ’ 


(2.11c) 


thereby requiring r < p if G" > 0 — a well-known requirement of thermodynamics — 
given that p > 0, r > 0, and 0 < a < 1. 

A Cole-Cole [17] diagram — a plot of G'(u ) versus G"(u) — is a very sensitive way to 
view anomalous relaxation phenomena. Figure 2.7 presents the normalized Cole-Cole 
diagram for the standard FOV fluid, with the a = 1 curve designating the response of 
a Maxwell fluid. Figure 2.8 presents a normalized Cole-Cole diagram for the standard 
FOV solid, which in this case translates the storage modulus of the FOV fluid by an 
amount ( t /p) q == 0.2, with the a = 1 curve designating the response of a Kelvin solid. 
The influence that the fractal order of evolution has on material response is readily 
apparent in a Cole-Cole diagram. Fractal order affects the extent of dissipation. 

Cole-Cole-type relaxations are naturally produced by random-walk models done 
on fractal lattices brought about by studying the motion of a particle in restricted 
geometries [37]. They also result from random- walk models done in fractal time on 
regular lattices, where the probability distribution is now a decaying power-law in 
time instead of the more common decaying exponential [43]. Random walks are used 
to establish the mean-square end-to-end distance of polymer chains in the various 
statistical theories of polymer physics, both for fluids and solids (e.g., cf. Douglas 
[25]). 


NASA/TM— 2002-21 1914 


26 







Chapter 3 

Continuum Mechanics: 
body-tensor fields 


In this report we use body-tensor fields (as defined by Lodge [68, 69, 70]) for deriving 
constitutive formulae, and we use Cartesian space-tensor fields (as commonly used 
throughout the literature) to solve boundary-value problems. In this regard we are 
free to select the tensor analysis scheme best suited for a particular task at hand, 
since one can readily map body-tensor fields into space-tensor fields, and vice versa. 
In this chapter we present the basic fields used in body-tensor analysis. In the next 
chapter we map these fields into Cartesian space, thereby producing tensor fields that 
are likely to be more familiar to the reader. 

Consider a continuum consisting of an infinite set of point particles, {fp}, also 
referred to as mass elements, that fills a domain, B, in 3-space (B C R 3 )- We call 
this set the body B. In any admissible body-coordinate system, B, defined over B, 
each particle <P in B is assigned a unique set of body coordinates, { = (IMM 3 ), 
E R, that are independent of time (i.e., B: fp — > cf- Lodge [69]). In this sense, 

body-coordinate systems have been construed as being convected coordinate systems. 


3.1 Metric Fields 

The distance dS separating any two neighboring particles — say, *P and <P' — in B can 
be quantified by using the covariant body-metric tensor, 7, of Lodge [68, 70] that he 
assigned to a Riemannian manifold with geometric measurement 

(dS 0 ) 2 = d£ ■ 70 • d£ and (dS) 2 = d£- 7- d£, (3.1) 

where 7(^3; t) is symmetric positive-definite and, most notibly, a function of time t 
with 70 := 7(^3 ;t 0 ), wherein t 0 is usually taken to be zero (0). This tensor field has 

a matrix representation of ^ (= [Yrcl; r,c= 1,2,3) in the coordinate system B with 
components y = y }i = y S calar dS(^]t) is the infinitesimal length-of-line of 

the contravariant vector d£(fp) := ^P t P / ; with dSo '■= dS^] t Q ) designating its gauge 
length. 
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A single dot placed between two tensor fields represents a contraction over one 
array index. Likewise, a double dot (colon) will represent a contraction over two array 
indices. 


3.1.1 Dual 

Because the metric tensor 7 is symmetric postitive-definite, it has a dual — the con- 
travariant (inverse) body-metric tensor, 7 _1 (^3;t) — whose matrix representation g -1 
(= [y rc ]]) has components = y- 7 ’ = y u (3,; t) in B such that 7 _1 • 7 = 5, where 
<5(fp) is the mixed idem tensor with a matrix representation of 6 (= [6£|) possessing 
components 6) that equal 1 whenever i = j, or that equal 0 whenever i / j, in every 
body-coordinate system. 

This dual body- metric tensor, 7 -1 , has two possible geometric interpretations. 
In the first interpretation, provided by Lodge [68, pg. 318], the distance 
measuring the height-of-separation between any two, neighboring, material surfaces — 
say, <r(£) = C and <r(£ + d£) = C + dC, where C and dC are constants — belonging 
to the same one-parameter family of non-intersecting surfaces, < 7 , in B is quantified 
via the contravariant body-metric tensor, 7 -1 , according to 


dC 

dH 0 


2 


da da 

dl'2° 'dl 


and 



da , da 

0 £'= 'd£’ 


dC = 



(3.2) 


where ( da/d £)(ip) is a covariant vector, independent of time t, that is normal to 
this material suraface, with contravariant vector d£ originating on one surface and 
terminating on the other. Scalar dH 0 := dH(i P; t 0 ) is the gauge length for this height- 
of-separation, while 7 q X := 7 _ 1 (^P; to)- 

A second geometric interpretation for the dual metric, provided by Truesdell [108] 
in an analysis done using general space-tensor fields, has a description of 

^ j~ l ^ 7 -! _ 

(dA 0 ) 2 = (d£ x d|) • • (d£ x dQ and (dA) 2 = (d£ x d|) • • (d£ x d$ , 

det q aet *y 

(3.3) 

where det 7 denotes the determinant of 7 , which is a scalar field with a tensorial weight 
of two ( 2 ), and therefore the areal metric tensor, (det 7 ) 7 -1 , is a relative field of like 
weight. Scalar cL4(*P;f) is an infinitesimal area-of-surface, with dA 0 := cL4(<P;f 0 ) 
designating its gauge area. This material surface contains neighboring particles fp, 
fP' and *P". The normal to this surface lies in the direction of a covariant vector field 
given by the cross product d£ x d£, which has a tensorial weight of minus one (— 1 ), 

wherein d£(f p) := fp*p' and d£(^3) := . 

The unit normal iy(^P; t) to element dA is given by udA := -v/detT-y d£x d£. This is 
an ‘absolute’ (without tensorial weight) covariant vector that, in body tensor analysis, 
is a function of time t through the presence of 7 (^ 8 ; t), because ||iy || 2 := ty- 7 -1 iy = 1 . 
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3.1.2 Rates 


The various descriptions for strain rate that are found in the literature can all be 
expressed in terms of derivatives of the metric tensor, 7, its inverse, 7 -1 , and its 
determinant, det 7. The time rate-of-change of the body-metric tensor, £>7, is defined 
classically through the limiting process 


7(^; t) - t(^P; f) 

DjO 3;f) = lim= , 


(3.4) 


while the time rate-of-change of the determinant of the body metric, D(det 7)(?P; t), 
is given by 

D (det 7) = (det 7) (tr £>7) , (3.5) 

wherein the trace, tr £>7, is computed as 7 -1 : £>7, consistent with the precepts of 
general tensor analysis. 

The first-order body-metric rates Dj($$\t) and £>7 _1 ((P;f) are not independent. 
Instead, they are related through the expression 

D7 -1 = — 7 _1 - (£>7) '7 -1 . (3.6a) 


Likewise, the second-order rates D 2 j[^\t) and D 2 7 x (<P;t) are related through 


£> 2 7~ x = 7 x - (2 (£>7) - 7 x - (£>2) - (-D 2 !)) ‘ 2 1 - ( 3 - 6b ) 

These identities are easily derived by applying the Leibniz product rule for differen- 
tiation to the expression 7“ x -7 = S, noting that D§_ = Q. Higher-order relationships 
can be acquired in like manner, but they are not needed in this work. 


Fractional order 


From the definition of Caputo differentiation (1.8a), one can compute fractal rates 
for metric evolution via the formulae* 

‘Because (“) / 0 whenever k > a given that k € N and a € K+ with a $ N, the Leibniz product 
rule (1.7) applicible to the Caputo derivative (1.3) leads to a more complex identity for -D? 7 ~ x than 
otherwise exists in the integer case (3.6a); specifically, for 0 < a < 1, 



m 7 


-1 


{to - t)-° 

r(l-a) 


(7 


1 

0 + 


00 / \ 

2->-E © 

*= 1 v 7 


(J fc -“7 _1 ) ' (£> fc 7) ' 7" 


where 7 “+ := 7 _ 1 (T; t 0 +). Consequently, there is no direct relationship between £>*7 and D *J. X - 
These rate fields are independent of one another. This result follows from Z>*( 7 _1 - 7 ) = = Q 
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where the fractal order of differentiation is restricted to the range 0 < a < 1. Un- 
like integer-order derivatives of the body metric (i.e. , and Z>y _1 (fP; t)), 

which are one-state fields in time t, fractal-order derivatives of the body metrics (viz., 
to, t) and to, t)) depend on all states traversed along the time in- 

terval [t 0 ,t], where we shall typically take t 0 to be zero in this work. 


3.2 Strain Fields 

Strain tensors are the preferred measures for describing deformations in solids, because 
solids have a quantifyible reference state. On the other hand, metric tensors are the 
preferred measures for describing deformations in fluids, becuase fluids have no unique 
state of reference. 

Strains are two-state fields that ideally possess four characteristic properties. The 
first property is: strain is a relative measure of deformation in that it vanishes when- 
ever its two dependent states are coincident. The second property is: strain is addi- 
tive and anti-symmetric in its dependency upon state. The third property is: strain 
exhibits tension/compression asymmetry; for example, axial extensions of A and A -1 
correspond to strains that are equal in magnitude but opposite in sign. And the fourth 
property is: strain is an absolute field (i.e., without tensorial weight) although, for 
the most part, this is really a requirement of convience. The second and third criteria 
actually quite restrictive. Only Hencky strain is known to satisfy all four of these 
criteria. 

The classic strain measures are defined below. The first two are tensor fields that 
relate to two distinct changes in length-of-line; the first tensor relates to a separation 
between neighboring particles, while the second tensor relates to a separation between 
neighboring material surfaces. The third strain measure is a scalar field that relates 
to the volume change of a mass element. 

3.2.1 Covariant 

The metric geometry of (3.1) can be rearranged in such a manner that (cf. Lodge [68, 
pp. 24-26]) 

{dS ) 2 - (dS 0 ) 2 = 2 ■ g ■ d£, g := ( 3 - 8 ) 

where e(fP; t 0 , t) is an absolute, symmetric, covariant, strain tensor. It has properties: 
i) tensor e vanishes in the reference state, e(^3;to,to) — 0; «) it is additive and anti- 
symmetric in its time agruments, e(f}}; to, t ) = e(^3; t 0 , t') + e(fP; t', t ) for all t' E [to, t], 
regardless of the extent of deformation; and in) it has no tensorial weight; however, it 
does not possess tension/compression asymmetry. Typically one sets to to zero, but 
for the time being, we shall leave it as to for clarity of discussion. 

The factor of 2 that appears in (3.8) is for historical reasons. Specifically, engi- 
neering strain is given by {l— lo)/l$ for the infinitesimal extension of a rod with length 
£(t) whose gauge length is t 0 := £(t 0 ). A normalized representation of the left-hand 
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side in (3.8) produces the relative strain measure 

Q dS f - (dS 0 ) 2 _ {dS + dS 0 ) (dS - dSo ) 

2(dSo) 2 2 dSo dSo (3 9 ) 

dS-dS 0 u 

« — whenever do ~ aoo, 

aoo 

demonstrating a consistency with the classic definition for engineering strain under 
conditions of infinitesimal deformation. 

Rates 

The covariant strain-rate tensor, is a one-state tensor field given by 

Dg=\Dj^ (3.10) 


because .D70 = 0. 

The covariant, fractal, strain-rate tensor, D“e(fP; to, t), depends on the path of 
straining incurred over the interval [to,t] of integration, and is given by 

(3-11) 

because H*7o = 0 (the Caputo derivative of a constant is zero). 


3.2.2 Contravariant 


In similar fashion to (3.8), the metric geometry of (3.2) can be rearranged in such a 
manner that (cf. Lodge [68, pp. 26-32]) 


dC 

dH 0 



(dC\ 2 = 2—-C-— C=H 7- 1 

\dH) di = d€’ =' 2l 2° 



(3.12) 


where £(^3;f 0l i) is an absolute, symmetric, contravariant, strain tensor*, which pro- 
vides another acceptable representation for strain. Like e, the strain £ is: i) a relative 

tProm (3.3), one can likewise define an alternative, symmetric, contravariant, strain measure as 

(<L4 0 ) 2 - (cL4) 2 = 2(d| x d|) •£• (d£ x d|), £ := l^det^o)!^ 1 “ ( det l)T _1 )> 

which is a relative measure of deformation in that /3(?P;<o>io) — and it is also additive and 
anti-symmetric in its time agruments because £(^3; to, t) = d(Ti ^o, f) + £(T; t',t) for all t' £ [f 0 , t\, 
regardless of the magnitude of deformation; however, it has a tensorial weight of two (2). It can 
be converted into an absolute field (i.e., without tensorial weight), but the outcome will violate the 
second desirable property of a strain field (viz., additive and anti-symmetric in state dependence). 

Because tensorial weights must equal amongst all additive terms in a tensor equation, and because 
all of the other tensor fields that we happen to use in constitutive development are absolute fields, it 
is therefore not practical to use 0 as a strain measure, even though it has the desirable interpretation 
of relating to changes in area-of-surface. 
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measure of deformation in that £(ty;t 0 ,t 0 ) = 0; ii ) it is additive and anti-symmetric 
in its time agruments because £(<P;< 0 ,i) = C(fP;to,^) + C(fP;i',t) for all t' E [t<j,t], 
regardless of the extent of deformation, and in) it has no tensorial weight; however, 
it does not possess tension/compression asymmetry. 

Rates 

The contravariant strain-rate tensor, D£(^3;t), is a one-state tensor field given by 

D i = = i” 1 ' ( D i) ' l" 1 ’ (313) 


because Z>y g 1 = 0. 

The contravariant, fractal, strain-rate tensor, £>“£(fP; ^o, t), depends on the path 
of straining incurred over the interval [to, t ] of integration, and is given by 

(3.14) 

because D^q 1 = 0 (the Caputo derivative of a constant is zero). Unlike the integer- 
order strain rates De and DC, there is no identity relating the fractal strain rate D^g 
to 

3.2.3 Dilatation 

To acquire a volumetric strain measure that is additive and anti-symmetric in its 
time dependency and exhibits tension/compression asymmetry, too, requires taking 
a different tact. Using the conservation of mass as our guide (viz., integrating Ding — 
— | tr D'y) leads to Hencky’s [50] definition for dilatation, A( f , P; t 0 , t), which is a scalar 
field given by (cf. Oldroyd [83]) 

A - | ln^det( 7 o 1- t)) = ln(p 0 /p) =\n(dV/dV 0 ). (3.15) 

This is our third classic strain measure. Like the prior strain measures, dilatation: i) 
vanishes in the reference state, Zl(fP;f 0 ,<o) = 0 ; ii) it is additive and anti-symmetric 
in its time dependence, <d(fP;f 0 , t) — A(ty]t 0 ,t') + A(ty]t' ,t) for all t' E [to,t], in- 
dependent of the magnitude of deformation; and in) it is without tensorial weight. 
Unlike the two prior strain measures, it iv) possesses tension/compression symmetry 
in that ln(dV 1 /dV 2 ) = -\n(dV 2 /dV l ). 

Scalar £>(fP; t) denotes the density of mass element with g 0 : = £>(^P; *o) being its 
gauge density. Scalar dV (*P; t) is the infinitesimal volume of mass element at time 
t , with dV 0 := dV(fP;f 0 ) denoting its gauge volume. Because det ^ 1 has a tensorial 
weight of minus two (— 2 ), while det 7 has a weight of plus two ( 2 ), it follows that 
det( 7 g 1 • 7 ) has no tensorial weight, and therefore A is an absolute scalar field. 
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Rates 

Dilatation A is actually defined via its rate, DA(ty\t), according to the expression 

DA := | trl >7 = -Din g = DlndV, (3.16) 

arises from the conservation of mass. It has a fractal rate of D“Z\(i}3; to, t), 

D ° A = f < 317 > 

is precisely the definition of a Caputo derivative (1.8a), assuming 0 < a < 1. 

Stress Fields 

is a linear mapping function defined by 

d(f> = 7£- v dA , and is constrained so that n = 7r T , (3.18) 

7r(^3;f) is the contravariant body-stress tensor introduced by Lodge [68, 70], 
is taken to be symmetric (i.e. , 7r = 7r T where superscript ‘T’ denotes the 
transpose). Its matrix representation nr (= |:t rc ]]) in the body-coordinate system B 
has components k ij = n jt = t). The resulting contravariant vector d(p(^P;t) is 

a differential force of contact acting on a material surface of infinitesimal area dA 
belonging to the mass element ^3 whose unit normal is given by u at time t. 

The differential area dA is assumed to be small enough (on a macroscale) that 
the differential traction vector d<f>/dA (which ratios force to area) is independent of 
its area, yet it must be large enough (on the microscale) that the contact force d</> 
exerted on area dA represents a statistical average taken over numerous inter-atomic 
and/or -molecular forces that comprise the mass element, thus granting us with a 
perspective, albeit vague, as to the (physical) size of a (mathematical) point in a 
continuum. 

There are times when it is preferrable to decompose stress into a sum of hydrostatic 
and deviatoric (i.e., traceless) contributions, which can always be done. Here the 
contravariant, deviatoric, stress tensor, 7 t(^ 3; t), is defined as 

7ir := 21 Tp7 -1 , and is constrained so that tr 7r = 0, (3.19a) 

where the trace, tr 7r, is computed as 7ir : 7. Prom this expression follows the definition 
for hydrostatic pressure, p(^3; t), which is a scalar field given by 

p ■= — |tr7r, (3.19b) 

wherein tr 7r is determined as 7r : 7. 

There are other times when it is preferrable to express stress as a contravariant 
extra-stress tensor, n(i]3;f), defined by 

II := 7r + p 7~ 1 , (3.20a) 


which 

where 

which 

3.3 

Stress 


where 

which 
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with the scalar p(<P; t ) being a Lagrange multiplier that is subject to an isotropic 
constraint of material incompressiblity; namely, 

det 7 = det7oi or equivalently, trL >7 = 0, (3.20b) 

recalling that tr £>7 = 7 -1 : Dy. The extra stress is not, in general, a deviatoric 
tensor; that is, tr II = 3 (p — p) need not be zero. 


3.3.1 Rates 


Like the metric-rate tensor, it is a straightforward matter to establish the stress-rate 
tensor, D 7 T, via 


D7t(!P; t) = lim 


IlOP;*) -2LpP;f') 

t - V 


(3.21a) 


so that the rate of deviatoric stress, Dn, becomes 


= D?g+ (Dp) 7 l +p(Dy 1 ), 


(3.21b) 


which itself is not deviatoric because D(^_ : 7 ) = D(0) = 0 implies that 

tr _D|| = ( ) : 7 = — 7 ; : (D 7 ), (3.21c) 

and where hydrostatic pressure evolves according to 

Dp = — : 7 + ^ : (£> 7 )). (3.21d) 

To compute the fractal stress-rate tensor, D°y_, one must solve the integral equation 

D;m to, t) = (Dk) CP; f) it', (3.22) 

where 0 < a < 1, with like expressions applying for and D+p. 

Derivatives of the extra stress n are handled differently, because the isotropic 
constraint p 7 -1 is actually a Lagrange multiplier and therefore it can be pulled 
outside the derivative. 


NASA/TM— 2002-21 1914 


36 



Chapter 4 


Field Transfer: 

Cartesian space-tensor fields 


Body-tensor fields, space-tensor fields, Cartesian space-tensor fields, and the map- 
pings between them, have all been carefully documented by Lodge in [68, 69, 70]. 
In Chp. 3 we presented the basic fields of body-tensor analysis. In this chapter we 
present an overview of Lodge’s mappings from the body into Cartesian space. These 
results are stated without proof. An useful artifact of any such transfer of field (i.e., 
mapping from the body into Cartesian space) is that the resulting spatial fields are 
objective (viz., frame invariant). We also present a section containing new results for 
the field transfer of fractional-order derivatives and integrals. 

The operation of field transfer makes it very plain as to whether a particular spatial 
field is Eulerian or Lagrangian. This characteristic of space tensors is affiliated with 
the time of field transfer. Eulerian fields result from a transfer of field at current 
time t , with this mapping being denoted by: body field space field ; whereas, 
Lagrangian fields result from a transfer of field at some reference time — say, t 0 (which 

we arbitrarily take to be zero) — as denoted by: body field i=^> space field. Without a 
knowledge of these field transfers, it is often difficult to ascertain whether a particular 

space field is Eulerian or Lagrangian. Detailing the underlying mathematics of t=> 
and t=^> is beyond the scope of this report. The interested reader is referred to either 
of the two texts by Lodge [68, 70]. 


4.1 Kinematics 

In contrast with the constructs of the prior chapter, continuum B can also constitute 
an infinite set of point places, {X 0 }, occupying a connected region in space, §, at 
some arbitrary time f 0 denoting its reference state. Each place X 0 relates to a unique 
particle in B and is given a label of X , which corresponds to the spatial position 
of X 0 (and therefore of ^3) in this reference configuration. Given an admissible, 
rectangular-Cartesian, coordinate system, C, defined over S, each place Xo in § is 
thereby assigned a unique set of spatial coordinates, X = (Xi,X 2 ,X 3 ), X* £ K, such 
that C : Xq — ^ X. 
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Later, at current time t (t > t 0 ), continuum B coincides with another infinite set 
of point places, {X}, that now occupies a different region in space, S. Each place 
X relates to a unique particle fp in B and is given a label of x with coordinates 
x = (xi, x 2 , x 3 ), which corresponds to the spatial position of X (and therefore of ip) in 
this current configuration. Consequently, C : X — » x where C is the same coordinate 
system used in the mapping C : X 0 — » X. 

Particle ip moves through space S with a velocity, y(f), of 


v := 


dx 

W 


(4.1) 


whose matrix representation is given by v = [v r ], with components, v* = dxi/dt, that 
are quantified in the Cartesian coordinate system C. 

The fundamental hypothesis of Cartesian continuum mechanics is that the motion 
at any location in the body is assumed to be sufficiently smooth in the sense that 
both 


5x 


F-SX and 


dF 

6y= 1 =-6X = 
at 


dg 

dt 



L ■ dx 


(4.2) 


exist, where g(t 0 ,t) := dx/dX defines the deformation-gradient tensor, and where 
L(t) '■= dv/dx defines the velocity-gradient tensor, neither of which is symmetric. 
They have matrix representations of F = |<9x r /c>X c ]] and L = \dv r /dx c J in the coor- 
dinate system C. The deformation gradient F is positive definite because, from the 
conservation of mass, 

0 < ^ = det(£) < oo, (4.3) 


and consequently, F -1 (f 0 ,£) = dX/dx always exists. In contrast, L is not positive 
definite, and as such, Lr l does not exist, in general. A subtle yet important fact is 
that £ and anchor to different locations; F anchors to X, while F -1 anchors to x. 

Particle fp changes its motion through space § with an acceleration, a(t ), of 

dv 

+ (4.4) 

whose matrix representation is given by a = [a,.], with components aj in coordinate 
system C quantified through the chain rule by a, = dvi/ dt + (dvi/ dx k )(dxk / dt) , where 
the repeated index k is summed over in the usual way. 

Position, x, velocity, y, and acceleration, a, are vector fields that establish kine- 
matic attributes belonging to a point in space. The deformation gradient, F, and the 
velocity gradient, L, are tensor fields that establish additional kinematic attributes 
belonging to a point in a continuum. 
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4.2 Deformation Fields 


In an Eulerian transfer of the various body-metric tensors into Cartesian space, Lodge 
[68, pg. 320] has shown that 


t t _ , \ 

7 l =^ L 7o §. I 

T 1 *=*Z. 7c 1 <=^ I j 


(4.5) 


where / is the unit tensor, and 5(3£; to, t) := £■ F T is the symmetric, positive-definite, 
deformation tensor of Finger [32], today commonly referred to as the left, Cauchy- 
Green, deformation tensor. Its reciprocal field B~ 1 (X\t( j ,t) = F _T • F ‘~ 1 is the actual 
deformation tensor that was introduced by Cauchy [14, pp. 60-69]. 

A Lagrangian transfer of these same body-metric tensors produces 

7 _1 i=^> c -1 , 


to j 

7o t=7 / 


t to 

7 o ' 


(4.6) 


where C(X 0 ]t Q ,t) := £ T -£ is the symmetric, positive-definite, deformation tensor 
of Green [46], today commonly referred to as the right, Cauchy-Green, deformation 
tensor. The inverse of this metric, C -1 (3Co; *o> £), is computed as F -1 - F~ T . 

t to 

It is apparent from the above results that the mappings t=> and t=^> are many- 
to-one. This consequence arises from the fact that Cartesian vector and tensor fields 
do not distinguish between tensorial kind (i.e., between covariant and contravariant 
indices in their coordinate transformation laws, because the Jacobian is restricted to 
be orthogonal for Cartesian fields), and as such, there is a loss of this information 
during these mappings.* 


4.2.1 Duals 

Finger [32] introduced both dual-metric tensors (viz., B and C -1 ) into the literature. 
It is well known that the fundamental metric tensors of the Eulerian and Lagrangian 
frames (i.e., B 1 and C, respectively) measure change in an infinitesimal length-of-line 
according to 

{dSof = dx-§T 1 -dx and {dS) 2 = dX- £■ dX. (4.7) 

Less known, and proven by Truesdell [108], is that the normalized inverse metrics 
B / detS and C -1 / det C _1 provide a like geometric interpretation; specifically, they 


‘General space-tensor fields do distinguish between kind in their coordinate transformation laws, 
and as such, field transfer between body-tensor fields and general space-tensor fields have mappings 
that are one-to-one. General space-tensor fields are not introduced in this report. The interested 
reader is referred to any one of the many excellent texts on the subject (e.g., Sokolnikoff [106]). 
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measure change in an infinitesimal area-of-surface according to* 


B 


cr 


(cL4 0 ) 2 = ( dx x dx ) • ~ - ■ ( dx x dx) and ( dA ) 2 = ( dX x dX ) • t • (dX x dX) > 

Q.6L D 


or equivalently, according to 


det C 


dA 0 

~dA 


B 


— n 


det B 


■ n and 


dA 

dA 0 


= N ■ 


C 


-l 


det C 


’ — > 


(4.8a) 

(4.8b) 


where n and N are the Eulerian and Lagrangian unit-normal vectors to a material 
surface, which relate to one-another via the pull-back formula NdAo = dX x dX = 
(det F)~ 1 F t • (dx x dx) = (det Z) _1 F T - ndA. Truesdell [108] closes his little-known 
paper with the following insightful theorem. 

“The elements of area suffering extremal changes are normal to the prin- 
cipal directions of stretch, and the greatest (least) change of area occurs 
in the plane normal to the axis of least (greatest) stretch; in fact, if the 
principal stretches dS/dS 0 satisfy > A 2 > A 3 the corresponding ratios 
dA/dA 0 satisfy A2A3 < A1A3 < A^.” 


4.2.2 Rates 

In an Eulerian transfer of field, Lodge [68, pp. 321-327] also determined that the 
various metric-rates of the body map into Cartesian space as 

D - 7 2D, 

D-f 1 i=4 -2£, 

where 

2(X;t) := \{L + L T ) 

is the symmetric rate-of-deformation tensor. The resulting rates, expressed 
some arbitrary tensor J, are defined by 

b=LA^l + L-L and | := L~ L‘ L~ L‘ , ( 41 1) 

which denote the lower- and upper-convected derivatives, respectively, of Oldroyd [83]. 
They reduce to Lie derivatives taken with respect to velocity y whenever dj_/ dt — g. 
The common contributing term in these two formute, 

| := § f + (V^) v, (4.12) 

tWe arrived at (3.3) by field transfer. Specifically, we mapped the tensor relations derived by 
Truesdell [108] from general space into the body, which is a one-to-one operation. Then by executing 
another transfer of field, this time mapping the formulae in (3.3) from the body into Cartesian space, 
which is a many-to-one operation, we arrived at (4.8a). 


(4.10) 
below for 


2ATo = 0 1 =^ B 1 = 0 
£>7 Q 1 = 0 t=4 B = 0 


(4.9) 
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is called the material derivative of J, which has a matrix representation of J = 
\d~S rc /dt + (dJ rc /dXk)wkl i n the coordinate system C. Here the vector operator V 
denotes the spatial gradient 8/ dx. 

A V 

The formulae B = 0 and B _1 = 0 of (4.9) can be rewritten as quasi-linear evolution 
equations for the Finger, B, and Cauchy, R~ x , deformation tensors; specifically, 

i = ]± R + R U and i _1 = -£> • JT 1 - iT 1 • 2, (4-13) 

where the resulting rate, when expressed for an arbitrary tensor J, is defined by 

l-=l-Kl + lK, (4-14) 

and is called the corotational derivative, which was introduced by Zaremba [110] and 
is usually credited to Jaumanny wherein 

t) := \{L - ^ T ) (4.15) 

is the skew-symmetric vorticity tensor. The quasi-linear evolution equation for Finger 
deformation in (4.13) lies at the heart of Leonov’s [63] viscoelastic theory. 

o 7 

The corotational J and lower-convected J deriviatives of any tensor J are related 
via 

i = ( 4 . 16 a) 

o A 

Similarly, the corotational J and upper-convected J derivatives relate according to 

1 = l+Ql + 12, (4.16b) 

where D is the rate-of-deformation tensor. From these identities, one quickly arrives 

o o 

at the evolution equations (4.13) for B and 5 -1 , given the field transfer results of 

<49) - . " 

We point out that rates 5(£; t 0 , t ), 5(X; t 0 , t ), 5 _1 (X; t 0 , t ) and J^ x (X; to, t ) are all 
two-state tensor fields. 

In Eulerian transfers of field where time-based derivatives of absolute, symmetric, 
second-rank, body-tensor fields are being mapped into Catesian space, the partial 
derivative of a covariant body tensor maps into a lower-convected derivative, while 
the partial derivative of a contravariant body tensor maps into an upper-convected 
derivative. 

* Contrary to numerous citations sprinkled throughout the literature, Jaumann did not publish 
the corotational derivative until 1911 [55], where it is given for both vector and tensor fields. An 
explicit definition of the corotational derivative does not appear in his book [54] of 1905. Zaremba 
[110] was the first to publish the corotational derivative of a tensor field, which he did in 1903. 
Nowhere (known to us) did Zeremba publish the corotational derivative of a vector field. 
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In a Lagrangian transfer of field, the body-metric rates map into Cartesian space 
according to 


D 7 DQ= 2g T QQ, 

£> 7 _1 ^ DQQ 1 = -2 £ -1 • a • Z“ T , 


-D 70 = 0 DI_ = 0 

= t > 

^o 1 = £ *=^ D l — Q, 


(4.17) 


where DC(X 0 \t 0 ,t ) and DC 1 (X Q ]t 0 ,t) are two-state tensor rates, and by DQ 1 we 
mean D(C~ 1 ) and not (DC) -1 , the latter of which does not exist, in general. 


Fractional order 

Because the process of integration commutes with the operation of field transfer [70, 
pg. Ill], it therefore follows from the preceeding formulae that the Caputo derivative 
(1.3) of the metric tensor, D°j(ty-,t Q ,t), maps into Cartesian space in a Lagrangian 
frame as D“ C(Xq\ t 0 , t) according to 


1 a r(t') „ 


1 c 1 <nv . 

D a ' v = — - — / = — dt 

= r(i - a) J t0 (t - t') a dt' 

fg-c- 1 /' 1 

*= r (i-c)C ((-(■)» de 


(4.18a) 


r (l 


In like manner, the Caputo derivative of the dual metric, D° 7 1 { f $;to,t), maps as 
D^C~ 1 (X 0 ]t 0 ,t) according to 


D“7 -1 = 1 [ 

= r(l - a) 4 


1 W* 


to 


D a C ~ 1 — 

*= r(i 


(t - t') a dt' 

1 r* 1 d£ -1 (f 0 ,i') , 


(■ t - t') a dt' 


r(i 




(4.18b) 


where the order of differentiation has been restricted to the range 0 < a < 1. The de- 
formation tensors present in these formulae are anchored to state to. These derivatives 
have units in time t of order t~ a . 

Given that a <b, the notation F(a, b ) implies dx(b)/dx(a), while QT l {a, b) implies 
dx(a) / dx(b ) , thereby establishing notation for a generalized deformation-gradient ten- 
sor. Changing the order of the arguments in the deformation gradient is equivalent 
to taking its inverse. So we adopt a notation for the deformation gradient wherein 
the numeric value of the first argument is always less than that of the second, and 
we invert the field, as needed, to satisfy this requirement. Whenever a = to, the 
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shortened notation of F<> := £(to,t') is used, and whenever a = t 0 and b = t, the 
accepted notation of £ := £(f 0 , t) is used, otherwise £(a, b ) is used. 

Following in this notation, let J f(a, 6 ) := £(a, 6 ) • £ T (a, 6 ) define a generalized 
Finger deformation and let C(a, 6 ) := £ T (a, 6 ) • £(a, 6 ) define a generalized Green 
deformation. In their arguments, the first state, a, denotes the initial (or reference) 
state; whereas, the second state, b, denotes the final (or current) state. 

In an Eulerian transfer of field, the Caputo derivatives D ° 7 and map into 

Cartesian space according to 


1 r* 1 d x (?) , 

DZi = -7 r / 7 -= — dt' 

* 7 r(i — a) J to (t — t 1 


a[ 

2 D := 


) a dt' 

1 f* 1 

r(i - a) k ( t - t'Y 


1 dB-\t\t) 

—— at 

dt' 


(4.19a) 


and 


D ^- l = 


F(1 


1^1 


to 
at 

" 2 - := f(l 


t'Y 


d 2r\t') 

dt' 


dt' 


r(i 


2— f 1 d MA dt ' 
- a) J t0 (t - t')° dt' 


(4.19b) 


with D first appearing in a paper written by Drozdov [27] but in a different notation. 
Notice that the deformation tensors in the integrands are now anchored to the floating 
state t', which is the dummy variable of integration, instead of to the fixed inital state 
t 0 . Also notice that the derivatives in the integrands (e.g., dQ(t',t)/dt') are taken at 
the reference state of the field (i.e., t') instead of at its current state (viz. t). 

The rate-of-deformation tensor, D, has no distinction between tensor kind; how- 

at a\ 

ever, its fractional counterparts have covariant-like, D, and contravariant-like, IJ_, 
constituents. 

The above expressions for the fractional rate-of-deformation tensors can be changed 
to anchor to the reference state f 0 by an appropriate application of the chain rule for 

differentiation. In doing so, the lower-fractal rate-of-deformation tensor, Q(X',t 0 ,t), 
of order a, with 0 < a < 1, can be rewritten as 


a[ 

D 


7-T 


F(1 








(4.20a) 


where the deformation gradients F T and F 1 can be taken outside the integral 
because their state dependence corresponds with the limits of integration. Likewise, 
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«r 

the upper-fractal rate-of-deformation tensor, D(X',to,t), can be expressed as 


- = ( r (l-a)/ t (t 

= -\g - {D^gr 1 ) -£ t 


f') Q =t' 


f ;, 1 • p(t') • £;, x dt' • £ 


(4.20b) 


From these results it is apparent that the Lagrangian fields D*C and D+C^ 1 are 
the fundamental (most basic) measures for fractional deformation rates in a spatial 
formulation in that they are actual Caputo derivatives. 

The Lagrangian and Eulerian, fractal, deformation rates of (4.18 & 4.19) are 
compatible (in the limit as a goes to 1 from below) with their first-order counterparts, 
which are: DC = 2 £ T • D • £ and DC -1 = —2 £ _1 • D ■ £~ T . 


4.3 Field Transfer of Fractional Operators 


We are now in a position to map fractional-order derivatives and integrals of any 
absolute, symmetric, body-tensor field into Cartesian space in both the Eulerian and 
Lagrangian frames of reference. 

We begin by considering an arbitrary, symmetric, covariant, body tensor, fi, and 
an arbitrary, symmetric, contravariant, body tensor, 77 , whose mappings into Carte- 
sian space are known. For purposes of illustration, consider a transfer of covariant 
field that is given by 


HOP;*) 


f t==> M(X; t ), 


such that N = £ T • M ■ £, 


(4.21a) 


and consider a transfer of contravariant field that is given by 


r^(<P; t) 


*==> K{% 


such that H_ — F 1 • G ■ £ T , 


(4.21b) 


where M and G are some arbitrary (but known), symmetric, Eulerian, tensor fields, 
with N and H designating their respective, symmetric, Lagrangian, tensor fields. 
In these transformations of field, the deformation gradient £ serves as a Jacobian 
of transformation between the Cartesian frames that pulls a known Eulerian field 
backwards, out of the Eulerian frame and into the Lagrangian frame. 

Lodge [68, pp. 321-327] has shown that the time rate-of-change D of a symmetric 
covariant tensor n maps into Cartesian space as 





M, 

DK, 


such that DN = £ T • M ■ £, 


(4.22a) 
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while the time rate-of-change of a symmetric contravariant tensor 77 maps as 

{ t * 

*7^ - such that DH = F~ 1 G-F- T , (4.22b) 

& DZ, - - - - 

where M and G are the lower- and upper-convected derivatives defined in (4.11), while 
DN and DH are simple first-order derivatives, as their notation suggests. The pull- 
back formulae in (4.22) can also be derived by differentiating the pull-back formulae 
of (4.21), as one would expect. 

Assuming that the above information is known to us, it is a straightforward matter 
to construct the special derivatives of interest. 


4.3.1 Derivatives 

Here we prove that the Caputo derivative £>“ (1.8a) of a symmetric covariant tensor 
H maps into Cartesian space as 


D> such that = (2>*M) •£, (4.23a) 

while the Caputo derivative of a symmetric contravariant tensor 77 maps as 

D a r) S *-’ such that D°H = F" 1 • (S)“ r G) ■ F~ r . (4.23b) 

= D?E, 

The Eulerian tensors £>“ L M and 2)° r G are objective rates of order a defined in (4.26), 
where notation ‘° 1 ’ is affiliated with covariant-like fields, while notation < “ r ’ is affili- 
ated with contravariant-like fields. Unlike D“N and D+ti, which are actual Caputo 
derivatives, derivatives 2 and are not Caputo derivatives, in a strict sense 

of the definition, which is why they are given a different notation. For this reason, La- 
grangian Caputo derivatives are considered the more fundamental of these fractional 
rates. 

In the derivations that follow, which prove the above mappings, a is restricted to 
the range 0 < a < 1. 

Proof: The Caputo derivative of symmetric body-tensor fields map into Cartesian 
space in a Lagrangian frame in an intuitive way. For example, for a covariant tensor 
H, its fractional derivative maps as 


D> = 


1 r* 1 d ^) df/ 

r(i - a) J t0 ( t - t'Y dr 

1=7 D a N = x _ r_j— 

*= r(i — a) J to (t — t') a 


dK(t 0 ,t') 

dt' 


dt', 


(4.24a) 
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whereas, for a contravariant tensor rj, its fractional derivative maps as 


1 r* l drj(t') 
■ ~ «) Jto V ~ * 


r (i 


t') a dt' 


<0 


D arr _ 1 /“ 1 

- F(1 - a) J t0 (t - f') Q 


d£{t 0 ,t') 

dt' 


(4.24b) 


dt'. 


These mappings reproduce standard definitions for Caputo derivatives of spatial fields. 
What is noteworthy about these mappings is that the argument list belonging to the 
spatial fields being differentiated under the integral sign is ( t 0 ,t '), designating a fixed 
referenced state at t 0 with differentiation under the integrand occurring at the second 
state t' . 

The mappings of these same Caputo derivatives (i.e. , and D°rj) from the 
body into Cartesian space in an Eulerian frame is a more subtle process. Before 
one can proceed, we must know how body-tensor fields from a past state map into 
Cartesian space in the present state, which they do according to 


a(t') ^Er T (t',t) -m(0 -r 1 (*',*)] 

- } t 0 < t' < t, (4.25a) 

ri(t')^F(t',t)-G(t')-F T (t',t) J 


where 

/x(f') M(0 and ij(t') G(f'). (4.25b) 

These transfers of field are distinct from the Eulerian transfers listed in (4.21). Here 
the deformation gradient E(t', t) serves as a Jacobian of transformation that pushs- 
forward a field that was anchored at some prior time t' into a like field that is now 
anchored at present time t. 

From the above formulas, it readily follows that the Caputo derivative of a sym- 
metric covariant tensor fx maps into Cartesian space in the Eulerian frame as 


D> = 


r(i 


1 —f 

- «) J t0 


1 x d^{t') t 
(: t - t') a dt’ ^ 


®* L M : = fn 


= rjTVj) l |rri yr F = r ( t ’’ *> ' M«') • £" *) 

= \T(l-a)J i 0 (t-tT dt’ dt ) = 


r T - 


(4.26a) 
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while the Caputo derivative of a symmetric contravariant tensor 77 maps as 


1 1 %(*') , 

D *= = r(i - a) J t0 (t - t'Y dt> dt 

’ '•= r(lb) £ jrhy I (S'. f > • m • E V, *>) * 

- mb) £ (br s' - ' () • s< ( '> • £ T «', 0 <*' (4 . 26b) 

= F ■ ( — — — r [ 7 - 1 ~t — F7, 1 ■ G(f') • F7 t dt'^j - F T 

= vr(i - a) J t0 (t - r) a =t = K =* y = 

- \T(l-a)J t0 (t-tr dt' ) = 

= F • (D°H) -E T - 


What is noteworthy about these mappings is that the argument list for the spatial 
fields being differentiated under the integral sign is in the defining integrals, 

designating a floating reference state at t' . In the equalities that follow the defining 
equality, this reference state is changed from floating (i.e., t ') to fixed (viz., t 0 ).O 


4.3.2 Integrals 

With the above derivations in place, it is a straightforward matter to prove that 
the Riemann-Liouville integral J a (1.2) of a symmetric covariant tensor /£ maps into 
Cartesian space as 

J a u [ X T 3Ql - such that J a N = F T - (T l M) ■ F, (4.27a) 

= J a K, - ~ “ 

while the Riemann-Liouville integral of a symmetric contravariant tensor 77 maps as 



^ such that J a H = F" 1 • (3 ar G) • F" T . 
J a !L ~ ~ 


(4.27b) 


The Eulerian tensors 3 QL M and Z al Q are objective integrals of order a defined in 
(4.29). Unlike J a N and J a H, which are actual Riemann-Liouville integrals, integrals 
Z a[ M and 3 Q, G are not true Riemann-Liouville integrals, in a strict sense of the defi- 
nition, which is why they are given a different notation. For this reason, Lagrangian 
Riemann-Liouville integrals are considered the more fundamental of these integrals. 

Proof: Like the Caputo derivatives, a Lagrangian map of the Riemann-Liouville 
integral for some covariant tensor n is a straightforward process whose outcome takes 
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the form 


JQ = r(a) 




(4.28a) 


while for some contravariant tensor rj it becomes 


1 f* l 

^ ja = = W) L ( * - i')— s( *" ,f,) *'• 


(4.28b) 


These mappings reproduce standard definitions for Riemann-Liouville integrals of 
spatial fields, and are valid for all a > 0. 

Using the mappings of (4.25), it follows that the Riemann-Liouville integral of a 
symmetric covariant tensor /x maps into Cartesian space in the Eulerian frame as 

t (429a) 

| -^-(r k!l¥^ Fjm ^ dt ')- F =' 

= l~ T - ( ja K) -ET 1 , 

while the Riemann-Liouville integral of a symmetric contravariant tensor rj maps as 
If* 1 

J a r) = -f- / t r\(H) df 

= r (a)J t0 (f-i') 1 -^ 

( N (4.29b) 

~ =£ - (fi = 

= £-(J q £)-Z t , 


for all a > 0.D 


4.4 Strain Fields 


Strain is not a unique concept in finite-deformation analysis. Here we present two, 
Cartesian, strain fields that relate to changes in length-of-line. The first pertains to a 
separation between material points, while the second pertains to a separation between 
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material planes. A third measure of strain accounts for changes in the volume of a 
material point. These concepts are presented in both their Eulerian and Lagrangian 
constructs. 

As it turns out, the strain fields used in the constitutive theories of Chp. 5 are 
different from the classic strain fields that are discussed in this chapter. Be that as 
it may, the strain fields presented in this chapter are of historical significance, and 
seeing how they arrive from field transfer should aid the reader’s understanding of 
how the spatial strain fields of Chp. 5 are arrived at. 


4.4.1 Covariant-Like 

The covariant strain tensor of (3.8) maps into spatial strain fields that are well known. 
In particular, in an Eulerian transfer of field, 

g := _ To) *=> A '■= |(Z - i _1 )> (4.30a) 

where A(X]t 0 ,t) is the Almansi [2] strain tensor, while in a Lagrangian transfer of 
field, 

g := 5(7- 7o) t=^> E •= \{Q ~ Z) » (4.30b) 

where g(£ 0 ;t 0 ,t) is the popular Green [46] strain tensor. 

The Almansi and Green strains are symmetric fields that measure a change in the 
distance-of-separation between a pair of neighboring material points. 


Rates 


From (4.22a), the time rate-of-change of strain g maps into Cartesian space as 



4 = g, 

Dg = | D£ = E T g-E, 


(4.31) 


7 

where A is the lower-convected rate of Almansi strain, and DE is the Green strain 

V 

rate, both of which are well-known results. Formula A = Q, can be rewritten in terms 
of the corotational derivative via (4.16a) as 


4 = U-E-A-AQ, (4.32) 

thereby producing a quasi-linear equation for the evolution of A in terms of its 
Zaremba-Jaumann derivative. 


Fractional order: The fractal, covariant, strain-rate tensor, D+g, of (3.11) maps 

into Cartesian space according to (4.23a) as 





£? l A = ZZ, 

= 1 0 - 2(0 ■&<*'. 


Dig = | d;q = r T - 2 ■ £ 


(4.33) 
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where It)* 1 A is the lower-fractal, Almansi, strain-rate tensor, and D+E is the fractional, 
Green, strain-rate tensor, both of order a wherein a is restricted to the interval 
0 < a < 1 . 


4.4.2 Contravariant-Like 

The contravariant strain tensor of (3.12) maps into Cartesian space in an Eulerian 
frame as 

£ : = Hio 1 -i -1 ) 2 ( 4 - 34a ) 

where Z{X]t 0 ,t) is a strain tensor that was introduced by Signorini [104], and whose 
appearance in the literature is scarce. A Lagrangian transfer of this same body field 
yeilds 

£ : = Kio 1 - t 1 ) ^ l : =UL- r 1 )’ ( 4 - 34b ) 

where YX%o',to,t) is a strain tensor whose origin is unknown to us, and that we shall 
simply refer to as a Lagrangian strain tensor. 

The Signorini and Lagrangian strains are symmetric fields that measure a change 
in the height-of-separation between a pair of neighboring material planes that are 
non-intersecting. 

Rates 

From (4.22b), the time rate-of-change of strain £ maps into Cartesian space as 




A 

where Z is the upper-convected rate of Signorini strain, and DY is the Lagrangian 

A 

strain-rate. Because of (4.16b), formula Z — D can be rewritten in terms of the 
corotational derivative as 

g = £ + + (4.36) 

thereby producing a quasi-linear equation for the evolution of Z in terms of its coro- 
tational derivative. 



Fractional order: According to (4.23b), the fractal, contravariant, strain-rate ten- 
sor, D“C, of (3.14) maps into Cartesian space as 


S)fZ = D, 




-l 


to 


D“Z= 


c*r 

D 


r T 


(4.37) 


dt', 


where T)fZ is the upper-fractal strain-rate tensor of Signorini, and D+Y is the upper- 
eth, fractional, Lagrangian, strain-rate tensor, both of order a wherein 0 < a < 1. 
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4.4.3 Dilational 

Using the preceeding results of this chapter, it is also a straightforward matter to map 
the dilatation A(t 0 ,t ) of (3.15) from a body formulation into a spatial formulation, 
thereby producing 


A = In ^det (7 q 1 -7) ^ i=4> A = In ^ 

(det B) 1/2 ) = In (det Z) 

(4.38a) 

in an Eulerian setting, and 



A — ln^det(7o 1 - 7 ) ^ l= =^ A = ln^ 

(det C) j = In (det JF) 

(4.38b) 


in a Lagrangian setting. There is no distinction between the Eulerian and Lagrangian 
descriptions for dilatation, as expected, because dilatation is an absolute scalar field, 
and is therefore independent of the Jacobian of transformation between these two 
frames. 

Rates 

A consequence of the above transfers for strain rate is that the time rate-of-change of 
dilatation (3.16) maps from the body into Cartesian space in an Eulerian frame as 

DA = | tr D-y i==>- DA = trZ), (4.39a) 

and it maps into Cartesian space in a Lagrangian frame as 

DA = | tr Dj DA = | £ _1 : DQ = tr Q. (4.39b) 

Again, there is no distinction between configurations because the transferred field is 
an absolute scalar. 

Fractional order: The fractal rate of dilatation, D*A, can be computed as 

D“A = — -i — - [ - 1 - (tr D)(t')dt', (4.40) 

* T(l-a)J 0 (f-f') a = JK ’ v ' 

which follows straight away from (4.39), or equivalently, from a field transfer of (3.17). 
This equation applies to both spatial frames of reference. 

Conservation of mass 

From the results of (4.38 & 4.39) it follows that the conservation of mass, whose rate 
form is 

D\xiq = — |trfl 7 t==4- Dlnp = — tr£>, (4.41a) 

integrates to 

f = ^=det£ (4.41b) 

in accordance with (3.15). 
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4.5 Stress Fields 


In an Eulerian transfer of field, Lodge [68, pp. 327-328] has shown that the body-stress 
tensor tt_ maps into Cartesian space as 


7 r T, (4.42a) 

where T(3C; t) is the Cauchy [14, pp. 42-56] stress tensor. The equation for deviatoric 
stress (3.19) therefore maps to Eulerian space as 


7r := 7T -f p'Y 


- 1 i==^ T :=T + pl 


i t , _ 

p := — ^ tr^ t=> p := — | tr 


(4.42b) 


where T is the deviatoric stress associated with hydrostatic pressure p in Cartesian 
analysis. 

In a Lagrangian transfer of field, the body stress maps into Cartesian space as 

K^E*, ( 4 . 42 c ) 

where P(X 0 ;to,t) := ^F _1 - T ■ F~ T is commonly referred to as the second, Piola- 
Kirchhoff, stress tensor, named after the independent works of Piola [85] and Kirchhoff 
[58]. Tensor P*( ^o! ^o, t) := F" 1 • T • F~ T is a Lagrangian stress tensor. § Both fields 
are symmetric, two-state, stress tensors. The deviatoric stress and the hydrostatic 
pressure of (3.19) therefore map into the Lagrangian frame as 

fK^E-E+fpgr 1 and fp^Ufp:=~lE-£, (4.42d) 

both of which are widely used definitions. 


4.5.1 Conservation Laws 

The physical law that impacts stress is the conservation of momentum. 

Balance of linear momentum 

This law yeilds a vector formula, called the stress equations of motion, that in an 
Eulerian frame is given by 

g a = V • X. + Qf_i (4.43a) 

where V • T is the divergence of stress, with a matrix representation of [<9T rfc /<9x fc ]] in 
coordinate system C, and where vectors a and / denote acceleration and body force, 
respectively. 

5 In the literature, the notation P* is often used to denote the first, Piola-Kirchhoff, stress tensor, 
If E -1 ’ I, which is not symmetric. The stress tensor that we call £* is seldom used. For us, tensor 
P* = F -1 • T ■ F- t is a natural choice, since it arises directly from a field transfer of the body-stress 
tensor 7r. 
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In a Lagrangian setting, the balance of momentum can be expressed as 

Qo q = Div(E-g) + Qof, (4.43b) 

where the divergence of stress is now taken with respect to the reference frame in that 
Div/7 has a matrix representation in C of |5H r fc/9X*| for some tensor ZZ Vectors a 
and / are the same as those appearing in (4.43a). 


Balance of angular momentum 

In the absense of couple-stress effects, which we assume herein, a balance^ of angular 
momentum requires that Cauchy stress X be symmetric, and therefore, X, X, X and 
P* must all be symmetric, too. 


4.5.2 Rates 

According to (4.22b), an Eulerian transfer of the body stress-rate tensor into Cartesian 
space leads to 

Dtt t=$- T, (4.44a) 

where T is the upper-convected derivative of Cauchy stress. A Lagrangian transfer of 
this same body field produces 

D^_ JU Dg\ (4.44b) 

where DP * = F _1 • T ■ F~ T > recalling that X = dX/ dt + ( VF ) • y — X ■ X — X • ■ 


Fractional order 

By constraining the fractal order so that 0 < a < 1, an application of (4.25) leads to 
an Eulerian transfer of field that produces 


a 1 [* 1 

T(1 - a) J 0 (; t - f') Q 


1 ^dt' 

dv 


S)“ r I = *) • W ■ *) dt ' ( 445a) 

= F • ( — — — - f — i-r- F-, 1 ■ T(t') ■ F~, t • F T , 


where the upper-fractal stress-rate tensor, T>* l X, of order a first appeared in a paper 
written by Drozdov [27], but in different notation. 
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A Lagrangian transfer of the same field leads to 


_ i i d*V) 

T(l-a)J 0 dt' 

1 r 4 i dP*(t 0 ,t ') 

to *= r(l — a) J 0 dt' 

*=* \ , rt , 




; t , T dt', 


where D^P* — F 1 - D^T ■ F T , which is compatible with its first-order counterpart 

A 

DP* = F~ 1 T ■ F~ t in the limit as a goes to 1 from below. 
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Chapter 5 

Constitutive Theories 


Continuum theories are presented in this chapter that are suitable for elastic and 
viscoelastic materials. Both compressible and incompressible constructions are dis- 
cussed. Each theoretical structure is reduced to cases of material isotropy and trans- 
verse isotropy. The resulting theories are mapped into Cartesian space in the Eulerian 
and Lagrangian frames of reference for use when solving boundary-value problems. 
Tangent moduli are also derived so these theories can be implemented into finite 
elements. Material functions are assigned to these theories in Chps. 7-??, thereby 
producing material models that can then be compared against experimental data. 

5.1 Integrity Bases 

The theory of invariants [107] is well developed as it pertains to our needs. In the 
case of a single tensor — say, A — having a matrix representation A in some admissible 
coordinate system affiliated with the underlying manifold, there exists an integrity 
basis that is comprised of three, irreducible, moment invariants: 

trA, trA 2 , trA 3 . (5.1a) 

This set constitutes an admissible, isotropic, integrity basis for any single, symmetric, 
3x3 matrix A. One advantage of using an integrity basis is that scalar fields can 
replace tensor fields as arguments in state functions. Another advantage is that an 
integrity basis leads to tensorial structure in a potential-based theory without having 
to introduce any ad hoc assumptions. 

None of the above three invariants can be expressed solely in terms of the other 
two; they are orthogonal to one another. Even so, this basis is not unique. The 
Cayley-Hamilton theorem — A 3 — (trA) A 2 -1- 1 ((trA) 2 — trA 2 ) A — (det A) I = 0, with 
invariants: trA, |((trA) 2 — trA 2 ), and det A — permits an exchange of the cubic 
moment (i.e., trA 3 ) with the determinant (viz., detA = |{trA 3 — trA[(trA) 2 — 
3 trA 2 }) when assigning an integrity basis. Furthermore, if matrix A is non-singular 
(i.e., detA ^ 0), then another application of the Cayley-Hamilton theorem allows 
the quadratic moment, trA 2 , to be replaced by the moment of its inverse, trA -1 . 
It is accepted practice to make use of both of these exchanges when assigning an 
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integrity basis for isotropic elasticity; specifically, trA, trA -1 , and detA are often 
employed because they lead to practicle models, even though they do not consititute 
an irreducible set. 

For the case of two distinct matrices — say, A and B — then in addition to those 
invariants that involve a single matrix only, one must also consider the coupled traces 

trAB, trA 2 B, trAB 2 , trA 2 B 2 , (5.1b) 

making a total of ten invariants that now comprise this irreducible integrity basis. 

When preference is given to one coordinate line above all other coordinate lines, 
then the resulting integrity basis is said to be transverse isotropic. Let vector a be 
a unit vector that is tangent to this coordinate line. Because there is no directional 
preference along a coordinate line, the integrity basis must be even in a, and as 
such, vector a can be represented by the symmetric dyadic matrix B := a <g> a with 
components B = [a r a c J, wherein <g> denotes the outer product between two vectors. 
Because a is a unit vector, it follows that a <g> a — (a Cg> a ) 2 = (a <g> a ) 3 = ... , and it 
also follows that tr a <g> a = 1, tr A • (a <g> a) = a • A • a, tr A 2 • (a <g> a) = a • A 2 • a, etc. 
Consequently, of the ten invariants in (5.1a k 5.1b) only five are independent; they 
are: 

trA, trA 2 , trA 3 , a • A • a, a -A 2 - a, (5.1c) 

which define the irreducible, transversely isotropic, integrity basis for matrix A. 

The intellectual process of selecting an integrity basis will lead to a constitutive 
theory. This is not a straightforward process because, although the irreducible in- 
tegrity basis is unique, there are numerous, other, orthogonal, integrity bases that are 
also admissible and therefore may be considered. Mathematics alone is incapable of 
selecting an integrity basis; physics must also be addressed. Once a theory is in hand, 
the process of assigning/deriving a potential function appropriate for the selected in- 
tegrity basis will produce a constitutive model. Only when such a model is available 
can the challenging process of characterization and (hopefully) verification against 
experimental data begin to take place. In physics, one cannot prove a theory/model 
to be correct. One can only, perhaps, disprove it, or more importantly, establish its 
domain of applicability. 

5.2 Elasticity 

Here we consider a mass element of density p(^3; t ) in a state of stress 7 r(^3; t ) that 

is undergoing an infinitesimal change in shape of := 7 (*P; t + dt) — 7 (^ 3 ; t) 

over an increment in time of dt. This induces a differential change in the work done, 
dW (SP; t), on the mass element (including the energetic contribution due to its change 
in kinetic energy) that is quantified by the formula [70, pp. 194-195] 

dW = ^K-dj. (5.2) 

2 Q— = 

From thermostatics, given an arbitrary reversible process, 7 -> 7 + dj_, the incre- 
ment of work done on a mass element minus the change in its kinetic energy (herein 
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denoted as dW ) equates with a change in the internal energy of the mass element. 
For an isothermal reversible process, the increment dW equates with a change in the 
Helmholtz free-energy of the mass element. 

Equation (5.2) leads to a potential-based constitutive equation for elasticity, 


f ^ = 


cm 

dy 


N 


n=l 


m din 
dl n dy ’ 


(5.3) 


which implies that (5.2) is an expression of the chain rule. The potential func- 
tion representing elastic strain-energy density, W (f$; 0, t ) = 2U(/i, / 2 > • • • , Jjv : ^P; 0, t), 
has arguments that make up an integrity basis comprised of N independent invari- 
ants I n {y 0 ,y,Xi,X2, ■ ■ ■ ,Xm-2- ^;0,£), n = 1,2,..., AT, obtained from M sepa- 
rate tensor fields: the metric tensors 70 and 7 along with M — 2 material tensors 
X m , m = 1, 2, . . . , M — 2. The tensor gradient dW/dy must be symmetric (i.e. , 
dW/d$ — |[9QU/9Y rc + dW/dy cr j in every body-coordinate system B: 7 -> 35) in or- 
der for it to satisfy the conservation of angular momentum. All anisotropic attributes, 
when present, are manifest through the gradient d%B/dy, which has units of stress 
on mass density. Unlike classical elasticity, there need not be a fourth-rank modulus 
tensor that assigns anisotropic characteristics in the theory of finite-strain elasticity; 
instead, the integrity basis introduces them through its gradients. 

When elasticity is described by a potential function then the material model is 
called hyper-elastic, and although there is no mathematical proof for its existence, 
Leonov [66] has provided physical proof for its necessity (there is no know proof of 
sufficiency): non-potential, finite-strain, elastic (so called hypo-elastic) constitutive 
relations can be constructed that create energy from nothing (i.e., they can operate 
as perpetual motion machines); hyper-elastic constitutive relations cannot be con- 
structed to violate thermodynamics in this way. This justifies referring to (5.3) as 
the fundamental constitutive equation governing elasticity. 

Theromostatics places loose constraints on what are admissible functional forms 
for the potential function 2IJ. This is an important topic in constitutive modeling, 
but it lies outside the scope of this report. 


5.3 Viscoelasticity 

One can ‘fluidize’ the above, elastic, constitutive laws, moving them beyond the realm 
of reversible thermostatics and into the domain of irreversible thermodynamics, by 
introducing a memory function that bestows a loss of remembrance onto past states 
[18]. The contributions that can be recollected from each past state are then summed 
over a loading history. Stress and strain are taken to be causal functions of time 
during this integration, and as such, are zero valued for all times prior to time t — 0. 
However, loading histories may be discontinuous at time t = 0, as will be the case in 
creep and stress- relaxation experiments. 

The well-known constitutive equation for linear viscoelasticity, restricted to in- 
finitesimal strains and rotations, can be written as the following Volterra integral 
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equation [47] 


T ij(t) = G ijU (t) E U (0, 0 + ) + G ijU (t - t') dt'. 


After an integration by parts, this expression becomes 


Tp-(t) — G{jkt( 0) E^(0, t) — 



E^(0 ,t')dt’, 


or equivalently, using the additive and anti-symmetric property of infinitesimal strain 
(i.e., Eki{0,t) — Eki(0,t') + Ekt(f,t), t' G [0,i]), it can also be written as 


T ij(t) = G ijU {t) E w (0, t)+ f M ijkt (t - t') E kl (t\ t) dt ', 

Jo 


where T*j = Tjj are the symmetric components of stress, E k t = E tk are the symmetric 
components of infinitesimal strain, Gj j k t — Gjiki = G ijt k are the symmetric compo- 
nents of a relaxation modulus, and M iju(t — t') := dG iju{t — t') / dt' are the symmetric 
components of a memory modulus. Memory fades if |My*/(£ 2 )| < |M tJ /^(i 1 )| for all 
h > h > 0. In classic viscoelasticity, the fourth-rank material functions G ijki and 
Mi j k t account for material anisotropy, when present. 

The first of the three formulations listed above requires the strain to be continuous 
and differentiable over time. The second and third formulations are less restrictive in 
that they only require strain to be continuous over time. The last two formulations 
differ in the moduli of their elastic terms, and they also differ in the states that define 
strain in their viscoelastic (integral) terms. It is the third expression of these three 
equivalent expressions that we choose to analytically continue from the infinitesimal 
into the finite. 

In all three of these classic formulations, strain is the controlled variable to which 
stress responds. Because the theory is linear, it can also be written so that stress is the 
control variable to which strain responds. But in our end application (finite elements), 
displacements are assigned to which forces respond, which motivates selecting strain 
for the cause and stress for the effect. 

Using classical viscoelasticity as our guide, and adopting the hypothesis of Kaye 
[56] and Bernstein et al. [7] (which they applied to viscoelastic liquids) wherein strain 
is replaced by the gradient of strain energy, as in elasticity theory, we therefore con- 
sider a class of K-BKZ - like viscoelastic solids that obey the constitutive hypothesis 


, ^ ~ di n (0,t) t ^dWdI n (t‘,t) J\ 

S = E - jp- + J, - * > W n ip- dt ) • (5 ' 4) 


and as such, produce a work increment (5.2) of 




71 = 1 


dW 

+ L 




07(0 
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where the viscoelastic strain-energy density = 20 (/ 1; / 2 , , In '■ '15; t', t) 

with invariants /«(7ts7, Xi, X 2 , • • • Jm- 2 : ^5 n — 1,2 is constrained 

such that 52B(fP; t, i )/$7 — 0 (there can be no strain when the two states defining 
strain are one in the same), and where 7c : = 7(15; t ') with t' G [0, <]. 

The n th relaxation function, ©„(£ — £'), which is positive and dimensionless, along 
with its associated memory function, 9Jt n (t - t ') := <9®„(£ — t')/dt', which is positive, 
has units of reciprocal time, and is monotone decreasing, do not specify anisotropic 
characteristics, when present, which is different from classical viscoelasticity. This 
task is relegated to the gradient of strain energy, d%B/d~f, which has units of stress on 
mass density. The fact that the memory functions, 27t n , are positive and monotonic, 
and that the state of integration, £', replaces the initial state, t — 0, as the reference 
state in the invariants of the integrand, are both in agreement with the hypothesis of 
a fading memory [18, 19]. 

It is because the memory function convolves with a first-order derivative (i.e., 
strain dW/dj), whereas the relaxation function convolves with a second-order deriva- 
tive (viz., strain rate d 2 W/dt chy), that motivated our selection of the K-BKZ - like 
constitutive structure presented in (5.4) as the phenomenological foundatation for our 
theory of viscoelastic solids. 

This type of construction (based on an elastic strain-energy density) is particularly 
attractive for modeling solids that are predominantly elastic with some viscoelastic 
attributes, like the materials of interest to us (viz., elastomers and soft biological 
tissues) . 

Equation (5.4), which is applicable for viscoelastic solids, is slightly different from 
the construct that one would use for viscoelastic fluids [71], which is 

g = ~Vs Dg 1 + 20o XJ / ~ 0^ dI 0(j) dt ' J trjD l~ 1 = 0 > 

wherein r) s represents solvent viscosity. Fluids are usually considered to be incom- 
pressible, and as such, the extra stress II replaces the body stress |[ in its construction. 
Furthermore, the viscous response of the solvent present in a viscoelastic fluid replaces 
the elastic response that is present in a viscoelastic solid. Finally, the lower limit of 
integration is moved from zero to minus infinity in the fluid theory, because fluids do 
not have a unique reference state. 

In contrast with the classic formulae of viscoelasticity, the constitutive theories 
presented in this chapter place no restrictions on the extent of either strain or rotation. 

Thermodynamics places loose constraints on what are admissible functional forms 
for the strain-energy density 2U, relaxation ® n , and memory 9Jl n functions. Like 
elasticity, this is an important topic, but it lies outside the scope of this report. 

5.4 Tangent Operator 

Commercial finite-element packages often require a user to supply a tangent operator 
(or modulus) for user-defined constitutive equations so that its solver can construct 
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an optimal stiffness matrix. Even though finite-element programs do not employ body 
tensors, it is still advantageous for us to construct this fourth-rank tensor in the body, 
and to then map it into the required frame of Cartesian space for later use in finite 
element analysis. 

The tangent operator £(*P;0,f) is a fourth-rank contravariant tensor defined by 

<£ := 2 — ^ ~ , so that d(^ k) = £ : (| dj) — £ : dg, (5.5) 

where e is the strain tensor of (3.8). This operator has symmetries (in a body- 
coordinate system B) of C jki = £? ikl = <& jik = <t jM , because stress zr and strain g 
are symmetric fields. In component form, (5.5) reads as = 29(^ K^)/dy ki- 
For the elastic constitutive equation of (5.3), the tangent modulus becomes 


N 


£ = 4£o 


dm d 2 I n 


N 




d 2 m dl„ 


n = 1 


dl n d'y dj ; dim dl n dj 

771=1 — 



(5.6a) 


which mandates an additional symmetry of = £ kh3 . The tensor outer product 
r) <S> rj has components in coodrinate system B. In most material models, 

including those of this report, d 2 m/dl m dl n = 0 for all m / n, thereby simplifying 
the above modulus so that it becomes 


N 


g = 4 £»o X! 

n= 1 


dm d 2 i n d 2 mdi n 

din dgd2 + dll dg 



(5.6b) 


In keeping with the assumption that there are no cross-coupling terms arising between 
the various invariants in the functions selected to represent strain-energy density, it 
follows that 


71=1 \ — — 


.dm d 2 i n 


dl n dj dy 


+ ^n(t) 


t')YYY^-dt' 

dl f 1 d 2 2L 

dll d ~1 JO 9I n 


d 2 m dl n dl n 


..d 2 m dl n ^ dl n 


(5.7) 


d'y d'y 


dt 1 


for the viscoelastic constitutive equation of (5.4). Some care is required here becuase 
In — 7 n (0, t) in the elastic terms; whereas, I n = /„(£', t) in viscoelastic (integral) 
terms. 


5.4.1 Stability 


For an elastic or viscoelastic solid to be stable requires that it be Hadamard stable, 
and as such, satisfies a condition of strong ellipticity where it becomes necessary and 
sufficient that 


€ 3kl x i y j x k yi> 0, V x,y. 


(5.8a) 
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A more useful condition (in accordance with 5.5), which is necessary but not sufficient, 
is that 

de:C:de = d{ s ^7£):de>0, V dg. (5.8b) 

This is a thermodynamic condition for stability, the so-called Drucker [28] stability 
postulate, which is not a physical law. 

The necessary and sufficient conditions for strong ellipticity are extremely complex 
and awkward to verify for any given model. It is for this reason that attempts have 
been made to obtain sufficient conditions for strong ellipticity that are of ample 
generality to be useful to the model developer. Of note are the works of Renardy 
[92], Leonov [65] and Kwon and Cho [61], all of which address viscoelastic liquids. 

Dissipative stability, pertinent for a stability analysis of viscoelastic liquids [65], 
has no role to play in the stability analysis of viscoelastic solids. 


5.5 Isotropic Elasticity 

An isotropic elastic solid can be formulated in terms of two, mixed, tensor fields 
when being described by body-tensor fields; in particular, consider A := and 
B := 35 -1 C 0 )* where one field is the inverse of the other. The conservation of mass, 
det( 7 _1 --y 0 ) — (p/p 0 ) 2 > 0, guarentees that 7 _1 -7o is non-singular, and consequently, 

invertible (and therefore, 7 q X - 7 is non-singular and invertible, too). As such, only 
three of the ten possible invariants in (5.1a & 5.1b) are independent. In what follows, 
we will postulate the existance of an integrity basis that is comprised of sums of like 
invariants taken from these two deformation tensors, and in doing so, we average their 
characteristics. 

From a mechanics perspective, there is no reason to choose the deformation vari- 
able 7 q 1- 7 over deformation variable 7 _1 - 7o, or vice versa. However, from a physics 
perspective, the molecular network theory of rubber elasticity (i.e., the neo-Hookean 
solid [53]) produces a constitutive equation described soley in terms of Jq ■ 7 [72]. 
This contrasts with the popular phenomenological model of Mooney [80] that uti- 
lizes both deformation fields in its description, and whose rigorous derivation from 
statistical mechanics continues to elude researchers. 

Rivlin [93, 94, 95, 96, 97, 98] and his students, Saunders [99] and Gent [39, 40, 41], 
were amoung the first to derive constitutive equations from an integrity basis for finite- 
strain isotropic elasticity, and to perform multi-dimensional experiments to seek out 
admissible functional forms for the strain-energy density. Reiner [91] was the first to 

*As is the norm in general tensor analysis, the trace of a tensor field (in this case, TJJ 1 and 
70) is achieved through a contraction with the metric 7 := y(<p,t), t > 0, or, when appropriate, 
its inverse 7" 1 := 7 _1 (ip ,t). It is because of this fact that gradients (i.e., dZB/dj) can arise from 
a strain-energy density W — 2U(ij;i = 1,2,..., IV) that may otherwise be considered as resulting 
from non-conservative sources, (e.g., one can introduce D7 as a state variable). This is one very 
important reason why we prefer using body-tensor fields for the purpose of deriving constitutive 
theories. 
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actually derive constitutive equations from an integrity basis. Today they are called 
Reiner fluids. 

Flory [34] proposed a multiplicative decomposition of the deformation field (defin- 
ing G := (det£) _1 / 3 F so that det G = 1) as a way to uncouple deviatoric and hydro- 
static responses, which has certain advantages from a computational perspective [105]. 
We choose not to make such a decomposition in this body of work because the in- 
variants we have chosen have gradients that are strain fields, not deformation fields 
as is typically the case, and therefore, it is not certain if the added complexity of 
introducing such a decomposition is warrented at this time. 

Adding respective invariants from the two deformation tensors, 7 q 1 -7 and 7 _1 -7o, 
each anchored to a different state, averages the effect of reference state that is tacitly 
held by these two measures for deformation and leads us to consider the following set 
of invariants as our integrity basis: 

h '■= \ (lo 1 : 2 + 1 : 2°) 

h : = w(io 1 -2'2o 1: 2 + 2~ 1 '2°'2 _1: 2 0 ) 
h := ^/det^o 1 - 2) + y^et^” 1 ■ ^o) 

whose gradients, * 

dh/dj = |(lo 1 -l- 1 -lo-r 1 ) 

dh/dx = I (lo 1 • 2 ' lo 1 - l -1 • 2° ' Z~ l ' 1 ° ' Zj 

dli/dg = \ (y/det^o 1 - 2 ) “ y / det(2“ 1 -2o)) 2 _1 


> , (5.9b) 


> , (5.9a) 


'Here we have used the following results to derive the gradients dl n / d'y and d 2 I n /d ■yd-f. 


d~f 


= 6 S £, with components in B of \ (6 \b\ + 5)5^) , 


_ 1 , 

d'Uj 



dg 



with components in B of 


3y kl 

dy ij 


= -Uv ik y ie + v u y jk )> 


and 


d \l\ 

d 2 


111 7 1 ' noti ng that llo 1 ’! 


Iio 1 ! Ill 


lil 

ii°r 
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are symmetric strain fields, and whose second derivatives, 


d 2 I x /d 7<?7 = |^7 1 El 7 1 - 7o- 7 1 + 7 1 • 7 o • 7 1 ® 7 X ) 
d 2 I 2 /dj^d g = | (^o 1 ® lo 1 + 1 _1 ' 1°‘ 1 _1 ® l -1 ' I 0- % X 

+ 7 -1 K1 7 _1 • 7o • 7 _1 • 70 • 7 _1 + 7 _1 ' 7o • 7 _1 - 7o * 7 _1 ® 7 _1 ) > , 

3 2 h/dgdx = I (\/ det (l _1 -lo) + y^det^o 1 - 7)) 7 _1 ® 7 _1 

+ l(\/ det ( 7 _1 'lo) - ^ det te 1- 2)) 7 _1 


(5.9c) 


appear in the tangent moduli of (5.6). At first glance, these invariants and their 
gradients appear to be somewhat complex but, as we shall soon discover, these def- 
initions produce relatively simple constitutive formulas in the Eulerian frame. The 
coefficients imposed on these invariants scale their associated strain measures so that 
they coincide with the definition of infinitesimal strain in the small-strain limit. The 
fact that a square root can be used directly in J 3 , but not in I x or I 2 , has to do with 
some unique properties of the determinant that are not shared by the traced 

Each strain measure in (5.9b) vanishes in the reference state; however, not one of 
these strains is additive and anti-symmetric in its dependence upon state; yet they 
all possess the desirable property of producing an asymmetric tension/compression 
response — a property that first appeared in finite-strain theory with Hencky’s [49] 
definition for strain, which is \ In ( 7 ^ • 7 ) = In A when expressed in terms of body- 
tensor fields [36]. 

Utilizing the gradients in (5.9b) obtained from the postulated invariants of (5.9a), 
the constitutive equation for an isotropic elastic solid derived from the work potential 

^Instead of using the mixed deformation tensors 7Q 1 - 7 and 7 _1 ‘7o as state variables, it may be 
preferable (especially for viscoelastic liquids) to use the mixed stretch tensors 

= := and = _1 = 

as state variables. From the identity A -1 • A = 5 , it follows that 


dA" 1 

#7 



and likewise, from the identity A • A = 7 0 : • 7, it follows that 


• 

d I 


A + A -5= = 70 13 d. 
— — 07 = - 


Like relations exist in the case of time derivatives. Only when a solution can be found for d£±/dj_, 
assuming that a solution does indeed exist, will it be possible to construct a theory using stretches 
instead of deformations as the state variables. 


NASA/TM— 2002-21 1914 


63 



in (5.2) has a stress response of 


f K = 2Qo ^233 1 1 (^o 1 - 2 ” 1 ' 2 ° ‘ l” 1 ) 

+ 2B i 2 |(2o 1 -2‘1o 1 ~2 _ 1 ‘2 0 ’2 _ 1 ’2°’1 _1 ) (5.10a) 

+ 2B 3 5 ('y/det^o 1 ' j) - ^det(2 _1 '2o))2 _1 )> 

whose trace defines a state of hydrostatic pressure (3.19b) 

fp = -|eo^,ii^o 1: 2-2 _1: 2o) 

+ 223,2 | (^o 1 • T ■ 20 1 : 2 ~ Qf 1 ' 2 ° ’ l 1 : 2°) (5.10b) 

+ 3223,3 |(y / det(2o 1 -2) “ ^/det^" 1 - ^o))) - 

that for materials whose bulk moduli are many times more stiff than their shear 
moduli becomes^ 

f -po22J, 3 (y / det(2o 1 -2) - y /det (2" 1 '2°))' (5.10c) 

Here 233, , := <9221(0, t)/dli, i = 1,2,3, are three material functions that can, at most, 
depend upon the three scalar invariants /i(0,f), I 2 ( 0, t) and / 3 (0,f). The invariants 
in ( 5 . 9 a) are defined as sums between like invariants taken from the two stretch 
tensors, 7 Q 1 - 7 and 7 -1 - 70 . Cureously, hydrostatic pressure p is described in terms 
of differences taken between these same, fundamental, stretch invariants. 

For the special case of an (ideally) incompressible material, the above theoretical 
structure reduces to 

g = 2<?o ^223 1 \ (^o 1 ~ 2~ l 'Z 0 ' l 1 ) + W ’ 2 8 (lo 1 - 1 ' lo 1 - 1" 1 ’ Jp- 2o- 1 _1 ) J , 

(5.10d) 

which is similar to setting 223, 3 = 0. Tensor g (= P 7 _1 + K) is the extra-stress of 
(3.20a), with p being a Lagrange multiplier that forces a constraint for incompress- 
ibility, det 7 = det 70 . 

The tensorial nature of (5.10) is fixed. Its structure is a direct consequence of 
the chosen integrity basis listed in (5.9). Only the three scalar functions 233, are 
adjustable. Their specification will turn this theory for elasticity into a material 
model, which is the primary topic of Chps. 7-??. 

§ It is for this reason that we do not find it necessary to decouple strain into separate hydrostatic 
and deviatoric parts. 
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5.5.1 Field transfer 

Using the results of Chp. 4, it is a straightforward process to map this theory for 
isotropic elasticity from the body manifold where it was derived into Cartesian space 
where it will be used to solve boundary-value problems. The main results from these 
mappings are outlined below. 


Lagrangian frame 

When the time of field transfer coincides with initial time t = 0, then the resulting 
Cartesian fields are referenced to the Lagrangian (or material) configuration. In this 
frame of reference, the isotropic elastic invariants of (5.9a) are computed as 


h — \ + tr Q 

/2 = ^((tr£ 2 ) + (trg- 2 )) 

h = Jdet C + yjdet CT 1 = det £ + det £ -1 


(5.11) 


where C (= £ T - £) is the symmetric deformation tensor of Green (4.6), often called 
the right, Cauchy-Green, deformation tensor, in which tensor £ (= dx/dX) denotes 
the deformation gradient defined in (4.2). 

An isotropic elastic solid whose invariants are so defined has a stress response 
(mapped from Eqn. 5.10a) of 

£ - 2p 0 (2D i \{L - Q- 2 ) + 2H.2 He - C- 3 ) + W, 3 \{detE~ det ET^Q- 1 ), (5.12a) 

with £ being the second, Piola-Kirchhoff, stress tensor defined in (4.42c). In the 
incompressible case, (5.10c) maps to space as the constitutive formula 

£* + pQ~ l = 200 (2U,1 1 (L - e -2 ) + 2B 2 l(£ - £~ 3 )) , (5.12b) 

subject to a constraint for incompressibility, det £ = 1. Tensor £* (= ^ £) represents 
the Lagrangian stress tensor of (4.42c) that, in the incompressible limit, is equivalent 
to the second Piola-Kirchhoff stress £ becuase q = f?o- 


Eulerian frame 

When the time of field transfer coincides with current time t, then the resulting fields 
in Cartesian space are referenced to the Eulerian (or spatial) configuration. In this 
frame of reference, the isotropic elastic invariants of (5.9a) are computed as 


£ = |(tr£ + tr £“ 1 ) 

*»=&((■ trS 2 ) + (trg- 2 )) 

/ 3 = yJdetB + yjdet £ -1 = det £ + det ET 1 


(5.13) 
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where g (= £ • g T ) is the symmetric deformation tensor of Finger (4.5), often called 
the left, Cauchy-Green, deformation tensor. Because the trace of a matrix product 
is independent of the order in which the product is taken, these Eulerian invariants 
equal their Lagrangian counterparts, which is why they are called invariants. 

An isotropic elastic solid whose invariants are defined by (5.13) has a stress re- 
sponse (mapped from Eqn. 5.10a) of 

f I = 2p 0 (2D lid- I _1 ) + au,2 I d 2 - i~ 2 ) + an 3 \ (det g - det ^- 1 )z) , (5.14a) 

with X. denoting the Cauchy stress defined in (4.42a). For the incompressible case 
(from a mapping of 5.10c), the stress response of (5.14a) reduces to an extra-stress 
response of 

z+pi= 2 6o (2u tl | d - r 1 ) + an, 2 id 2 - r 2 )) > ( 5 - i4b ) 

where the Lagrange multiplier p forces a constraint for incompressiblity, detF = 1. 
Special care must be taken when implementing this constitutive equation (in fact, 
when implementing any nearly incompressible material model) into finite elements 
so as to avoid possible ill-conditioning of the stiffness matrix leading to a potential 
locking of the mesh caused by an overconstrained displacement field. 

It is here, in the Eulerian frame, that the full simplicity of our constitutive theory 
for isotropic elasticity (derived from the integrity basis of 5.9a) is most apparent. 

Equation (5.12a) can be pushed forward into the Eulerian frame producing (5.14a) 
by contracting (5.12a) with F from the left and F r from the right. Conversely, (5.14a) 
can be pulled backward into the Lagrangian frame producing (5.12a) by contract- 
ing (5.14a) with £ _1 from the left and F~ T from the right. These are appropriate 
push-forward/pull-back mappings for contravariant-like tensor fields. Because these 
mappings apply to Z and g*, and also to ^ T and P, it is sufficient to obtain a single 
transfer of field from the body into Cartesian space (hereafter, we will provide map- 
pings into the Eulerian frame) from which the corresponding formulation in the other 
(viz., Lagrangian) frame can then be readily acquired by applying the appropriate 
(i.e., pull-back) mapping between these two Cartesian configurations. 

Generalized strain fields: The elastic theory presented above is constructed with 

generalized strain fields of the type 

I (n) : = ik(l" - I~")» " e R, (5.15a) 

with dilatation being approximated by the scalar measure 

e := |(det£ — det^ -1 ), (5.15b) 

where instances n = 1 and n = 2 are the strain fields present in (5.14) so that^ 

f Z = 2f? 0 (2D, 1 1 (1) + 2D 2 i (2) + 2D 3 eg), (5. 16a) 

^Another constitutive theory that may prove useful is 

f 1= 2p 0 (oI {1/2) +6| (1) +celj, 
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(5.16b) 


with a hydrostatic pressure of 

f p = -|<?o (SET i tr g 1] + 2D, 2 tr£< 2 > + 3 2B 3 e) « — 2^ 0 2U, 3 e, 

that for incompressible materials becomes 

z + P L = 2po (an, 1 ! (1) + 22J 2 £ (2) ) , (5. 16c) 

satisfying a constraint for incompressiblity, detF = 1. We will explore the capabilities 
of these elastic constitutive formulae in the remaining chapters of this report. The 
notation E (7^ D n E ) does not imply a derivative of order n, as it is often understood 
to mean. 

Equation (5.15) is distinct from the generalized strain fields introduced by Doyle 
and Ericksen [26], where they define E^ ee := C" — £). Equation (5.15) is similar 
to, but still distinct from, the generalized strain fields of Bazant [6], where he defines 
:= T_( C n — CT n ). Actually, E ^ = R ■ E$ ■ R T with g denoting the orthogonal 
rotation tensor (i.e. , R T ■ R = /) gotten from a polar decompostion of g such that 
F = V • R = R ■ U, where V and U are symmetric positive-definite tensors called 
the left- and right-stretch tensors, respectively. It is easily verified that Q = ¥? and 

Q = g 2 . 

Taylor expansions (expressed as series in infinitesimal strain) of any Doyle-Ericksen 
strain field and of Hencky [49] strain, \ lng = ln^7, are coincident only in their linear 
terms; they differ in all other terms in their expansions. In contrast, Taylor expan- 
sions of any Bazant strain field and of Hencky strain are coincident through their 
quadratic terms, differing thereafter [6]. Taylor expansions of (5.15) and | ln^ = In Y. 
are likewise coincident through their quadratic terms, independent of n. This makes 
E a reasonable approximation for logarithmic (true) strain, without having to deal 
with the complexities that otherwise accompany Hencky strain [51]. 

The strain fields of Hencky, Bazant, and Eqn. (5.15) all possess the desirable 
property of tension/compression asymmetry; for example, stretches of A and A” 1 in 
uniaxial extension have strains that are equal in magnitude but opposite in sign. 


Tangent operator: In an Eulerian frame where, for example, tangent moduli are 

constructed in updated-Lagrangian finite-element codes, the tangent operator of (5.5) 
maps into Cartesian space as^ 

(fl) A =£:a (5.17) 


where a, b and c are material functions of some set of invariants yet to be determined. Such a theory 
could be constructed from a work (or strain energy) potential if one knew how to solve 


dg • 

dC 


U+U 


9 g 

dg 


7 ISI /, 


which here is represented in terms of Lagrangian fields. 

I4n a Lagrangian frame, the tangent operator of (5.5) maps into Cartesian space a s dg = g: dg, 

where P is the second, Piola-Kirchhoff, stress tensor and g is the Green strain tensor. The Lagrangian 
and Eulerian tangent moduli have components that relate according to Cijkt — ^Tm^Jn^kp^JqCrnnpq- 
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where D is the rate-of-deformation tensor defined in (4.9), T is the upper-convected 
rate of Cauchy stress given in (4.44a), and where the fourth-rank tensor Q_ is the 

t 

Eulerian tangent operator obtained through the field transfer <g i=> Q_. 

The tangent operator for the isotropic elastic solid of (5.14) is given by 

g=4Q 0 (m,i\({l®g~ 1 ) + +W n {g {1) ®K {1) ) 

+ 237 2 1 ((i B g) + (if 1 B § _1 ) + (I H if 2 ) + Qf 2 H Z)) + 237, 22 (f (2) 0 1 (2) ) 

+ 2H, 3 e(L 12 Z) + (2U.3 | h + 227 33 e 2 ) {[ 0 Z)) , 

(5.18) 

where 227 j = dW/dR and <27 ^ — d 2 W/dIidIj, i,j = 1,2,3, and where tensor prod- 
ucts MHN and M<g>N have matrix representations of MEI N = ||(Mj fc N J 7 + lVMNfj fc )J 
and M 0 N = in the coordinate system C. This tangent modulus was 

obtained by substituting (5.9b & 5.9c) into (5.6b) and mapping the resultant into 
Cartesian space. Not all of the terms in (5.18) will be contributing in any given 
model; furthermore, additional terms will be required if the function representing 
strain-energy density has a cross-coupling of I m In, m ^ n, in its definition. 


5.6 Isotropic Viscoelasticity 

With a general structure for isotropic elasticity in place, it is a straightforward process 
to extend this structure to a theory for isotropic viscoelasticity using the constructs 
of (5.4), which is based upon the K-BKZ [7, 56] hypothesis. As in the case of isotropic 
elasticity, the stress response has three constituent parts, one related to each invariant, 
that combine to quantify the total state of stress according to 


f K = 2q 0 (©r(t) 237i(0,f)| • 2<r 1 _1 ) 

+ J STti (t - t ') 237, i {{, t ) \ - T 1 ’ Jt' ' 'l l ) dt ' 

+ © 2 (f) 227 2 (0, t ) | (^o 1 • 2 ’ lo 1 ~ 7T 1 ■ 1<> ■ T 1 • 2° ' 2” 1 ) 

+ J 97^2 (t- t ') 237, 2 {t', t) | • X I*' 1 ~ 1 1 ' 1* ' l 1 ' l 1 ' ' 2 _1 ) dt ‘ 

+ ® 3 (t) 237,3(0, t) l^-^/det^o 1 ' 7 ) ~ y'det^- 1 -^))^” 1 

+ J 27 l 3 (t -2) 237,3(2, t) det (g-, 1 • 2) - yj det (2 _1 ■ X ) ) dt ' > 

(5.19) 

where the invariants of (5.9a), adjusted for time interval [t’,t], lead to constitutive 
formulae that satisfy the required constraint of <9227(f, f)/07 = g whenever t' = t. 
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There are three, independent, scalar- valued, relaxation functions 0;, i = 1,2,3, 
in this consitutive equation whose associated memory functions are computed as 
— t') := d<£>i(t — t')/dt'. The three 0* are dimensionless, positive and bounded. 
Their memory functions - f) have dimensions of reciprocal time, they can 

possess a weak singularity at t' = t, they must be positive, and they must decrease 
monotonically with a sufficiently rapid rate of decay towards zero as t goes to infinity, 
all in accordance with the hypothesis of a fading memory [18, 19]. The fact that the 
memory function can be weakly singular at t' — t is offset by the required constraint 
that d%D(t' ,t)/dy — 0 whenever t' = t, so the integrand as a whole remains regular 
over the entire interval of integration. 

One possible simplifying assumption that can be considered is to set the relaxation 
moduli 0t and 0 2 equal to one another, indicating that shear flow is governed by a 
single relaxation mechanism, independent of the order of strain. 


5.6.1 Field transfer 

Using the results of Chp. 4 and, in particular, Eqn. (4.25), the stress response of the 
isotropic viscoelastic solid listed in (5.19) has an Eulerian description of 


f I = 2 eo (©1 (t) 20 1 (o, t) £ (1) (0, t) + jT (t - t') 20 1 (if, t) £ (1) (if, t ) dt' 

0 2 (t) 20 2 (0, t) £ {2) (0, t) + f m 2 (t - tf) 20 2 (t', t) t ) dt' (5.20a) 

Jo 

0 3 (t) 20,3(0, t) e(0, t) L + 2 Jl 3 (t - O 20, t) e(t', t ) dt' lj , 


+ 

+ 


that, for materials whose bulk modulus is much greater than its shear modulus, 
produces a hydrostatic pressure that (approximately) satisfies 


fp= -2p 0 (®3 (t) 20,3(0, t ) e(0, t) + jT W 3 (t - t') 20 3 (t', t ) e(t', t) dt'^j . 


For incompressible materials, this constitutive description simplifies to 


(5.20b) 


Z+pI = 2q 0 (ei(t) 20 1 (0, t) £ {l) (0, t ) + jT JOt! (t - t') 20 1 (if, t ) £ {1) (if, t) dt' 
+ 02 (t) 20 , 2 (0, t) | (2) (0, t) + Jm 2 (t- If) 20 , 2 (If, t ) g 2) {t', t ) dt'^j , 


(5.20c) 


subject to a constraint for incompressiblity: either det£ = 1 or tr£> = 0. Here 
E(t',t) := dx(t)/dx(t') = £(0, t) ■r 1 {0,t') and g{t',t) := £(?, t) ■ £ T {t', t) so that 
g n \t',t) = ^{W{t',t) -§r n {t',t)) and e{t',t) = |(det£(t',t) - detiT^t)). This 
is the theory for isotropic viscoelasticity that we employ for solving boundary-value 
problems. 
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Tangent operator 


The tangent operator (5.7) is composed of two parts for each invariant — an elastic 
part and a viscoelastic part — whose sum produces the overall tangent modulus for 
the isotropic viscoelastic solid of (5.20); specifically, 

C = C e + C v , (5.21a) 


such that upon substituting (5.9b & 5.9c) into (5.7), and then mapping the resultant 
into Cartesian space, one obtains an elastic part that when expressed in the Eulerian 
frame is given by 


(5.21b) 


g e = Ag 0 (^> 1 (t) ^20,i + (ZT^Z)) +Wn{E {1) 

+ © 2 (t)^2iT2|((ii3|) + (zr 1 ^- 1 ) + (z^jr 2 ) + 

+ 2 v 22 (g 2) ®g 2) )^J 

+ <S 3 (t) (an, 3 e(£ H z) + (211,3 \ h + 2U,33 e 2 ) (Z ® Z)) j , 
and a corresponding viscoelastic part that is given by 

g v = *qo(J o 2j i,{t - + (r^t'.t) mL))dt' 

Jo 

+ (ZH|- 2 (t',t)) + (§r 2 (t',t)®l)^dt' (5.21c) 

+ [ m 2 (t-t')<m i22 {t , ,t)(g 2) (t',t)®g 2 \t',t))dt' 

Jo 

+ J m 3 {t-t')W, 3 (t',t)(e(l®l) + \l 3 (t',t)(l®L))dt' 

+ j m 3 {t - t') 2 V 33 {t', t ) e 2 (t', t)(L ®ljdt'^j , 

with additional terms required whenever the function representing strain-energy den- 
sity has a cross-coupled dependence between its invariants. 


5.7 Transversely Isotropic Elasticity 

We now address materials composed of a single family of fibers, or showing a single 
preferred orientation. Let the contravariant vector a o := a(f)3, 0) be a tangent to this 
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material line of anisotropy at a particle ip in the reference state of time t — 0, and 
let this vector have unit length in this state so that it satisfies 


«o-7o- «o = 1 


(5.22a) 


At current time t, the fibers along this material line will have stretched by an amount 
A according to 

A 2 = c*o - 7 • c*o- (5.22b) 

Any material obeying such a symmetry property is called transversely isotropic, be- 
cause the transverse material plane (defined by a normal vector that is coaxial with 
the covariant vector 7 0 • c*o) remains isotropic. Virtually all biological tissues are at 
least transversely isotropic. Many are orthotropic, and therefore exhibit two prefer- 
ential directions, but we will not address this or any other higher level of anisotropy 
in this report. 

In order to construct an elastic theory for transverse isotropy, in accordance with 
(5.1c), we need to introduce two additional invariants to the three isotropic invariants 
that already exist and which are defined in (5.9). The two additional invariants that 
we propose to use are defined by 


h := |ao- + 2°' 2 1- 2°) 0 

h := (7 ' lo 1 ’ 1 + 2° ' 2 _1 ‘ 2° ‘ Z ~ l ' 2°) 


(5.23a) 


whose gradients are anisotropic strain measures given by 

dh/dj^ = |((a<} ® « 0 ) -2 -1 ‘2 °' '2° 'Of 1 ) 

dh/d^= ^((a 0 <g>« 0 ) -Tl? _ 2 _1 '2 0 ' (^ o<8> ^°) ■2 0 '2 _1 '2°'2 _1 ” ’ 

+ To*" 7 ' (a 0 <8>a 0 ) - 7 _1 ‘ 7o- 7 -1 ' 7o- (ao®a 0 ) -7o-7 _1 ) ^ 

(5.23b) 


and whose second-order derivatives, present in the tangent moduli, are 

d 2 h/dj : d2 = i( 2 _1 b 2 _1 '^ 0 ' «o) ‘ 2 0 ' 2 _1 

+ 7 _1 ■ 70 • (a o 0Q O ) • 7o - 7 _1 Sy; -1 ) 
d 2 h/djdj= j^((a 0 ® a 0 ) ^ lo 1 + 70* ® («o®ao) 

+ 7 _1 E3 7 _1 - 7o- (ao <8> a 0 ) • 70 ■ 7 -1 ' 7o* 7 -1 
+ 7 -1, 7o- (ao®a 0 ) -7o-7 _1 ® 7 -1 ’ 7o- 7 _1 ' 

+ 2~ 1 '2 0 ' (^0 <s» oc 0 ) •2 0 ’2" 1 '2 0 '2 _1|a 2" 1 

+ 7 _1 Kl7 _1 -7o-7 _1 -7o • (a 0 ® a 0 ) -7o-7 _1 
+ 7 _1 ' 7o • 7 _1 7 _1 • 70 • (a 0 ® «o) • 7o • 7 _1 

+ 7 _1 '7o- 7 _1 '7o • (ao®ao) '7o-7 _1 ^7 _1 ) _ 


(5.23c) 
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Collectively, the five invariants of (5.9a & 5.23a) constitute an admissible integrity 
basis for transversely isotropic materials.** Gradients dli/d 7 and dh/d 7 are strains 
of first order; whereas, gradients dl^/d^ and dl^/d^ are strains of second order. 

From these five invariants, one arrives at the constitutive equation for an elastic 
solid with transverse isotropy derived from the work potential in (5.3); it has a stress 
response of 


f K = 2 ^o^ 2 U,i *) 

+ 22J 2 1 (^o 1 • g ■ go 1 - t 1 • 2° ' 2 1 ’ 1° ' 2f 0 

+ 2n 3 |(^/det(2o 1 -l) - ^det(2 _1 -2o))l -1 

+ 20,4 |((«o 0 «o) -7 _1 '7o- (a 0 ® ao) 

2XJ 5 ^((a 0 ®a 0 ) - 7 - 70 1 - 7 _1 -lo- («o®ao) • ^o- 1”*' 2°' 2 

+ 7 q 1 • 7 • («o®ao) — 7 -1 ' 7<r 7 -1 ' 7o - (a 0 <8>«o) ■ 2?' 2 *) 


-1 


(5.24) 


There is an apparent symmetry in tensorial structure between the contravariant fields: 
70 1 and a 0 ® a 0 , and between the covariant fields: 70 and 70 • (a 0 ® Qc 0 ) ■ 7 o- 


5.7.1 Field transfer 

Before one can proceed with a mapping of this constitutive equation from the body 
into Cartesian space, it is necessary to determine how the unit vector a 0 of (5.22) 
maps into space. In a Lagrangian transfer of field, let the Lagrangian unit vector a 0 
be defined by the field transfer 


a 0 t=b a 0 

ao- 7<r ao = 1 £o - ao = 1 >, 

, o ^0 ~ ■» 2 

QLo' — ^ 0<r C • a 0 = A 


(5.25a) 


**It was our desire to construct meaningful, anisotropic, finite-strain measures that led us to 
choose invariants that are themselves a sum of like invariants constructed from fields with opposing 
variance (e.g., one constituent invariant arises from the contravariant form of a state variable, while 
the other arises from the covariant form for that same state variable). Constructing meaningful, 
isotropic, strain measures is not a difficult task; there are many admissible choices. However, there 
are few admissible choices to select from when attempting to construct meaningful, anisotropic, 
strain measures. Equation (5.23b) presents one such pair of acceptable strain fields. 
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so that in an Eulerian transfer of field one gets 

t . 

a.o >=>■ Xa 

a o • 7o • ao = 1 Q. * |[ _1 ' « = A -2 > , 
«o*7‘C*o — A a • a = 1 


(5.25b) 


and therefore a — X 1 E ■ ao, or equivalently, ao — XF 1 • a. Material curves defined 
by the trajectory of a 0 are indicative of fiber direction. After a deformation £, these 
fibers have a new direction of a, and they have stretched by a factor of A. Vectors a 0 
and a are the same vector fields that were utilized by Spencer [107, pg. 13]. 

To the invariants of (5.11 or 5.13) we add (gotten by a field transfer of Eqn. 5.23a) 

I 4 = \ < 1 q- {C + C~ l ) ■ ao and I 5 = {Q 2 + Q 2 ) • Ho> (5.26) 

which are different in form from the invariants used by Spencer: floT'^o and flo'£ 2 '^o- 
Affiliated with the invariants of (5.26) are the two, anisotropic, Eulerian, strain fields 


4 (1) := 1 A 2 ( (a <0 a) - £ _1 ■ (a <g> a) ■ 

'■= jo X 2 (^((a< 8 > a) ■ g + g- (a® a)) > , 

- BT 1 ■ {(a 0 a) ■ £ _1 + ST 1 ■ (« 0 «)) • JT 1 ) 

which result from a transfer of the fields in (5.23b). Recalling that 

X 2 =g.Q-C-a 0 and a = X~ 1 F-a 0 , 


(5.27a) 


(5.27b) 


allows these anisotropic strain fields to be recast as 
A (1) - ao) ® (E-a 0 ) - {ET T -a 0 ) ® (Z~ T '^o)) 

4 (2) = A ((£ ■ a 0 ) <0 (g ■ E ■ ao) + {E • E ■ «o) ® {E ■ «o) 

- (r T -«o) ® {sr 1 ■ e~ t ■ ao) - (i' i -£ _T -«o) 0 (r T -«o)) 


(5.27c) 


These are still Eulerian strain fields, it is just that they are represented in terms of 
the material vector ao instead of its spatial variant a. In the deformed state, vector 
F ■ a 0 is tangent to the material line of anisotropy, while vector g~ T - «o is normal to 
the plane of isotropy, as illustrated in Fig. 5.1, and these vectors need not be coaxial. 

A compressible, transversely isotropic, elastic solid whose invariants are so defined 
has a stress response (mapped from Eqn. 5.24) of 

f Z = 2 qo (2Dr 1 E {1) + 2B , 2 E {2) + 233 3 e L + 233 4 4 (1) + 22J 5 A {2] ) , (5.28a) 

that in the incompressible case becomes 

I +pL = 2 Q o (an, 1 1 (1) + an, 2 1 (2) + 2n 4 4 (1) + 2 r, 5 4 (2) ) , (5.28b) 

where the Lagrange multiplier p forces a constraint for incompressiblity, det £= 1. 
These are the formulations for transverse-isotropic elasticity that we investigate in 
the remaining chapters of this report. 
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Tangent operator 

To the isotropic, elastic, tangent operator listed in (5.18), one must add an anisotropic, 
elastic, tangent operator (obtained from substituting Eqns. 5.23b & 5.23c into Eqn. 5.6, 
and then mapping the resultant into Cartesian space), which is given by 


g a = 4<?o 22J 4 A 2 (z El Qf 1 • (a 0 a) ■ JT 1 ) + (JT 1 • (a 0 a) ■ g _1 ) El lj 

+ ^ 22J 5 A 2 El (a 0 a) + (a 0 a) El B 
+ /El • ((fl®a) • B~ x + 5 _1 - (a 0 a)) • 

+ ((a 0 a) • g- 1 +§T 1 ‘ (a0 a)) El/ 

+ 2T 1 El (if 1 - (a 0 a) - I" 1 ) + (if 1 - (a 0 a) - I" 1 ) EIZT 1 ) 

+ 233,44 (4 (1) 0 4 (1) ) + 233,55 (4 {2) 0 4 (2) ) . 


(5.29) 


Again, additional components will need to be included into this expression if there 
exists a cross-product dependence between the invariants in the function representing 
strain-energy density that is being considered. 
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5.8 Transversely Isotropic Viscoelasticity 

In accordance with our elastic theory for tansverse isotropy (5.24), an application 
of the K-BKZ - like constitutive hypothesis for viscoelasticity (5.4) produces a vis- 
coelastic constitutive relation governing the stress response of a transversely isotropic 
material that has five constituent parts, one associated with each invariant, and is 
given by 

f K = 2*?o 2B !(0, f) | (^0 1 - Z 1 - 2?' I 1 ) 

+ J* - 1') \ (77 1 - T l ■ jt' • 2 _1 ) dt ' 

+ 02 (t) 22T 2 (0, t) | ( 27 1 • 2 • lo 1 - 1 _1 • 2° ■ l -1 ■ 1° ' l -1 ) 

+ J m(t - 1 1 ) 2 t) | (77 1 • 2 ■ zi - 1 ~ z 1 • Z* ■ Z ' 1 ■ Z}' • l " 1 ) dt ' 

+ 0 3 (i)2B, 3(0,t) |(^/det(27 1 ‘2) ~ ^/det ^” 1 - 2 0 ))2~’ 

+ J m 3 {t - t')W 3 (t\t) ■ 1) ~ ^/ det ( 2 “ 1 - 2 t '))^ / z " 1 

+ 04(020,4(0,0 I ((a 0 < 8 >a 0 ) - 7 _ 1 - 2 0 ' (^o®«o) - 2 °' 2 _1 ) 

+ J 3Jt 4 (t-0®,4(<' ) Oi((a('®«(')-2 1 2<' - (at- <g>atO 

+ 0 5 (O 2 U, 5(0,0 s((«o®«o) ■ 2 ' 2 o % ~ Z 1 ' Zs>' («o® So) ■ 2 o , Z~ 1 'Z?'Z~ 

+ To 1 ‘ Z ' («o® «o) -7 _1 -7o-7 _1 '7o- (ao®«o) ' 2°' Zj 

+ [ m 5 (t-t')W 5 {t' } t) 

Jo 

X Ye^ioLt' ® a f ) • 7 • 77 1 - 7 _1 • 7f • («<' ® «t') * 2«' ' Qf 1 ' ZJ ’ ' Of* 

+ 77 1 • 7 • (a* ® ««0 - z 1 ' Zy ' 7 _1 • 7 t* • (« t' ® «<') ■ 1 *' ‘ 2 ” 1 ) ^ ) ’ 

(5.30) 

which contains isotropic viscoelasticity (5.19) as a special case. 

A possible simplifying assumption that one may choose to consider would be to 
set 0i = 02 and 0 4 = 0 5 , thereby implying that relaxation behavior does not 
discriminate between the first- and second-order strain fields, yet it does discern a 
difference in the relaxation behaviors along the strong direction and in the plane of 
isotropy. 

Another simplifying assumption to consider, which is applicable for at least some 
biological tissues [101], is that 0i = 0 4 and 0 2 = 05, with the implication being 
that relaxation behavior does not vary between the strong direction and the plane 
of isotropy, even though the elastic response can possess a strong anisotropy; yet, 
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unless ®i = 0 2 , too, (and therefore 0 X = 0 2 = 04 = 0s) the material can discern 
a difference between the relaxation behaviors due to first- and second-order strain 
effects. 

A more useful assumption would be to eliminate the higher-order invariants / 2 and 
h that introduce nonlinear capabilities, provided that experiments, or the boundary- 
value problems to be solved, justify making such a simplification. 

5.8.1 Field transfer 

The vector of anisotropy in the reference state of integration, a t i , maps into Cartesian 
space in this floating frame of reference as follows, 


t' 


QLt' a t ' 

t' 

QLt' ■ It' ■ QLt' — 1 *=> Lit' ■ a? = 1 
QLt ' ' 7 ‘ QLt' = A 2 (f', t ) 1 =^- a? ■ Q(t', t ) • a v — A 2 (t', t) 

and consequently, it maps into the Eulerian frame as 

OLt> t==> A (t 1 , t ) a 

QLt' ■ It' ■ QLt' = 1 a ■ ,t)- a = A 

QLt' • 7 • QLt' — A 2 (t', t) !==>■ a ■ a = 1 


(5.31a) 


(5.31b) 


where A(f',f) denotes the stretch along the fiber axis over the time interval 
Furthermore, the anisotropic strain fields are given by 

A {1) (t',t) = \ A 2 (f',f)((a<® a) - (§T\t',t) • (a 0 a) -g" 1 ^,*))) 

^ A 2 (f',f)^(((a®a) -§{t',t)) + (B(t',t) ■ (fl®fl))j ; 

-JT' 1 (f',t) ■ |((a®a) + (|| _1 (t',i) • (a® a))) •£ _1 (* , 1 *) 


wherein 


with 


A — a t '- C(t',t) ■ a t i and a = A 1 {t' ,t) F(t' ,t) • a t ’, 


A 2 (0, f') = a 0 • C(0, t') ■ a 0 yielding a t ' — A x (0, t') F(0, i') • a 0 , 
which generalize the formulas of (5.27). 


(5.32a) 

(5.32b) 

(5.32c) 
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The stress response of the transversely isotropic viscoelastic solid given in (5.30) 
has an Eulerian description of 


Z = 2 g 0 (t) 2 B i(0, t ) K {1) ( 0 , t) + ^ m (t - t') 20, i (f , t ) £ (1) (t', t) dt' 

+ 0 2 (t) 2 B 2 (o, i) i (2) ( 0 , f) + r an 2 (t - f') 20 2 (t', i) i (2) (*', t) dt' 

Jo 

+ 0 3 (t)20 3 (O,t)e(O,i)Z + [ W 3 (t - t')W 3 (t' ,t) e{t' ,t) dt' L (5-33a) 

Jo 

0 4 (f)2B 4 (O,t)4 (1) (O,t)+ f 2Jl 4 (t - t') W 4 (t', t) 4 (1) (t', t) dt' 

Jo 

0 5 (t) 20 5 (o, t) 4 ( 2 ) (o, t) + J* m(t - 1') 20, 5 (*', t) 4 ( 2 ) (f', i) dt' J , 


+ 

+ 


that, for an incompressible material, reduces to 

l + p L = 2 qo (® x (t) an, i (o, *) i (1) (o, t) + jf * ^ (t - t') 20 x (*' , t) £ (1) (<', t) dt' 

+ 0 2 (t) 2U )2 (o, t) i (2) (o, t) + f an 2 (t - t') 20 2 (t', t) z {2) (*', t) dt' 

Jo 

+ 0 4 (t)2n 4 (o,t)4( 1) (o,t)+ [ m 4 {t-t')m 4 (t',t)A (1 \t',t)dt’ 

Jo 

+ 0 5 (t) 2 XJ 5 (o , f) 4 (2) ( 0 , t) + r m 5 {t - 1 ') 2 Dr 5 (^, f) 4 (2) (t', t) dt ') , 

do / . 

(5.33b) 

subject to a constraint for incompressiblity: either det|: = 1 or tvQ = 0. These 
are the formulations for transverse isotropic viscoelasticity that we exercise in the 
remaining chapters of this report. 


Tangent operator 

To the two parts already listed for the isotropic, viscoelastic, tangent operator in 
(5.21), two additional parts must be considered to account for transverse isotropy, 
such that 


c = c e + c ea +c v + c va , 


(5.34a) 
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where the isotropic components Q_ e and Q are found in (5.21b & 5.21c), and where 
the anisotropic components are given by 

g ea = 4p 0 C 0 4 (*) \\ 2 (l® (if 1 -{a® a)- if 1 ) + (if 1 • (a® a) • if 1 ) 0 [j 


+ 2U 


',44 



+ ® 5 (t) 2X1 5 ^ A 2 ( B Kl (a ® a) + (a ® a) IE B 


+ g- 1 B (if 1 ■ (a ® a) ■ if 1 ) + (if 1 -{a® a)- g~ l ) 0 g~ l 
+ i El ^if 1 • ((a ® a) ■ BT 1 + Br x ■ ( a ® a )) • B -1 ^ 

+ (fi -1 • {{a ® a) • if 1 + if 1 • (n ® a)) • if 1 ) 0 A 

+ 0 5 (t)2H,5 5 (4 (2) O4 (2) )Y 


and 


(5.34b) 


g ua = 4p 0 ^ 2Jt 4 (t - t') 253 4 (f i) | A 2 (f f) (i B (if Y, <) • (a ® a) • g~ l (t', t )) 

+ QfY,*) • (a <8) a) -ifft', f) Bi)dt' 

Jo 

+ ^ 97i 5 (t - t') *) A Y'> 0 (lY f B (a 8 a) + (a ® o) B |(f t) 

+ if V, f) S (g~\t', t)- (a® a)- g~\t\ t)) 

+ ( ; 1 {t',t) ■ {a®a)-g~ l (t',t)) B|fY,t) 

+ ZK I (jfY, t)- ((a® a) ■§r 1 {t',t)+g- 1 (t',t)- (a® a)) * ^ _1 (t / , t)) 

+ (jfff f) • {{a®a)-g- 1 (t , ,t)+g~ 1 {t',t) ■ (a® a)) -JfY, f)) B^dt' 

J*m 5 (t - t')W 55 {t\t){A {2 \t\t)®A [2 \t\t))dA, 


+ 


(5.34c) 

with additional terms required if the function representing strain-energy density hap- 
pens to depend on cross products between its invariants. 
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Chapter 6 

Finite-Strain Experiments 


Although many characterization experiments are ID in concept, they are 3D in reality 
and often require 3D analysis, even for the simplest of deformation histories, espe- 
cially when finite strains are being considered. We restrict our discussion to those 
experiments whose stress and strain fields remain spatially uniform (homogeneous) 
throughout the specimen. Furthermore, deformations are assumed to be controlled to 
which forces are measured as reactions, while temperature is held constant. This re- 
striction of controlling deformations and measuring forces, although imposed herein, 
is by no means a requirement of characterization experiments. Shear-free extensions 
and simple shearing are the deformations considered. 

Deformation rates and strain rates are known because they are controlled. Force 
rates and stress rates, on the other hand, are usually not measured. Whenever stress 
rates are present in a constitutive equation, be they of integer or fractional order, 
a system of differential equations ensues that needs to be integrated to arrive at a 
predicted stress that can then be contrasted with the observed stress. This requires 
that stress rates be handled different than strain rates for purposes of characterization. 


6.1 Shear-Free Extensions 

A deformation is said to be shear free if there exists a convected coordinate system 
(e.g., B : ip -» Jj) that is always orthogonal [70, pg. 81]. In this class of experiments, 
a uniform material element in the shape of a unit cube is deformed into a rectangular 
parallelepiped causing tractions to set up on the various material faces, dependent 
upon boundary conditions, as illustrated in Fig. 6.1. These experiments are techno- 
logically important for solids because of their relative ease of execution. They are, 
however, extremely difficult to execute on fluids, where they have only met with some 
success when testing polymer fluids of high molecular weight. 

6.1.1 Kinematics 

The deformation just described locates an arbitrary particle in test sample B with a 
set of spatial coordinates X — say, X ~ (Xi,X 2 , X 3 ) — in some rectangular-Cartesian 
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( 0 , 1 , 0 ) to 



Figure 6 . 1 : Shear-free extension of a material element. 


coordinate system C at a reference state of to- Later, this same particle is located at 
a different place with coordinates x — say, x = (a^, X2, x 3 ) — in the same coordinate 
system C, but now at the current state of t. These two coordinate pairs relate to one 
another through the motion 

£1 — A x X\, X2 — A2 X2 and X3 — X3X3, ( 6 - 1 ) 


where A* is the principal stretch in the i th direction. 

The deformation gradient tensors, F and F _1 , therefore have components of 



fAi 0 0] 


X 1 0 0 ' 

F = 

O 

V 

to 

0 

and F 1 = 

0 A2 1 0 


1 

CO 

'C 

O 

O 
1 


0 0 A3 1 


whose polar decompostion has a rotation matrix of 


10 0 
0 1 0 
0 0 1 


leading to stretch tensors V and U with identical components of 



Ai 

0 

o' 


Ai 

0 

o' 

V = 

0 

A2 

0 

and U = 

0 

A2 

0 


_0 

0 

A 3 _ 


_0 

0 

A 3 _ 


whose inverses are trivially 



'Af 1 0 0 ' 


rv 0 0 1 

V" 1 = 

0 A2 1 0 

and IT 1 = 

0 

tH 

1 CN 

O 


0 0 A3- 1 


L 0 0 A 3 


(6.2) 


( 6 . 3 ) 


( 6 . 4 ) 


( 6 . 5 ) 
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Obviously, there is insufficient information in these components to exploit the full 
richness of these kinematic fields. 

The velocity gradient, D and the rate-of-deformation, £, tensors have like com- 
ponents of 

D In Ai 0 0 1 plnAi 0 0 

L = 0 Din A 2 0 and D= 0 DlnA 2 0 . ( 6 . 6 ) 

0 0 D In A 3 J 0 0 D In A 3 _ 


There is no vorticity (i.e., W = 0) because the eigenvectors of strain are fixed in B. 

<*1 . 

The lower-fractal rate-of-deformation tensor, D, has a matrix representation of 


ai 0 0 
D = 0 IXfD^Xl 0 

0 0 |A 3 - 2 D“AI. 


(6.7) 


af 

while the upper-fractal rate-of-deformation tensor, D, has a representation of 


t r-|A 2 D“Ar 2 0 0 

D = 0 -|A 2 D?A 2 2 0 

0 0 -PIoja , 2 


( 6 . 8 ) 


a) of 

such that in the limit, as a goes to one from below, both Q and U. yield Q as expected. 


6.1.2 Deformation Fields 

In the Eulerian frame, the contravariant-like, Finger, deformation tensor, and the 
covariant-like, Cauchy, deformation tensor, £ -1 , have components of 

A 2 0 Ol \Xf 0 0 ' 

B = 0 A| 0 and B _1 = 0 X^ 2 0 . (6.9) 

0 0 A|J 1_ 0 0 A3 2 . 

In the Lagrangian frame, the covariant-like, Green, deformation tensor, Q, and its 
rate, DC, have matrix representations of 

‘A 2 0 Ol pAxDAi 0 0 

C = 0 A 2 0 and DC = 0 2A 2 DA 2 0 , (6.10) 

0 0 A 3 0 0 2 A 3 DA 3 . 

while the contravariant-like inverse, C _1 , and its rate, DQC 1 , have representations of 

A " 2 0 0 1 T- 2 Aj; 3 DAi 0 0 

c~ x = 0 A 2 2 0 and DC -1 — 0 - 2 A 2 3 DA 2 0 

0 0 A 3 2 0 0 — 2 A 3 3 DA 3 

( 6 . 11 ) 
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The Caputo derivatives of these deformation fields have components of 



'D? A 2 

0 

0 ] 

’D« A," 2 

0 

0 

d:c = 

0 

D? A 2 

0 and D^C- 1 = 

0 

D? A 2 2 

0 


0 

0 

D* A 2 J 

0 

0 

^a 3 - 2 _ 


(6.12) 

for 0 < a < 1. 

The fact that the deformation fields B and C have identical components is a 
direct consequence of the deformation gradient F being diagonal. In general, the 
Finger and Green deformation tensors have matrix representations that are different 
from one another, but not for shear-free extensions. 

6.1.3 Strain Fields 

With the deformation tensors and their rates now known, it is a simple matter to 
quantify any strain or strain-rate field of interest. In what follows, we determine the 
five Eulerian strain fields used in our models. 

The three isotropic invariants defined in (5.13) have values of 

h — \ (A? + A 2 + A 3 + 2 + A 2 2 + A 3 2 ) 'j 

h = ^ (^1 + A 4 + A 3 + Aj 4 + A 2 4 + A 3 4 ) > , (6.13) 

h = AiA 2 A 3 + A^A^Ag - 1 J 

which are observed to be insensitive to tensions versus compressions, an achievement 
that other classic invariants can only obtain through squaring, whereby making them 
even functions. Associated with these invariants are the isotropic strains from (5.15) 
that in component form are given by 

„ riM- A r 2 ) 0 0 1 

EW = I(B-B- 1 )= 0 |(A|-A 2 - 2 ) 0 , (6.14) 

0 0 i( A 3-A3 2 )_ 

"§(A 4 — A^ 4 ) 0 0 

J3 (2) = |(B 2 — B~ 2 ) — 0 !(A^-A 2 - 4 ) 0 , (6.15) 

. 0 0 §(a 4 -a 3 - 4 )_ 

and 

e = |(|F| - |F| _1 ) = i(A,A 2 A, - Ar'AJ-Aj 1 ), (6,16) 

where the mass densities ratio as 

= A 1 A 2 A 3 , (6-17) 

in accordance with the conservation of mass. 

Figure 6.2 plots stretch against strain for the deformation of simple extension 
where curves for E^, E^ and (infinitesimal) engineering strain are presented. What 
is observed is that E^ is a very good approximation for engineering strain to quite 
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Strain Measures in Uniaxial Extension 



Figure 6.2: A comparison of various strain measures for the ID extension of a rod. 


large extensions, while is true nonlinear measure of strain. All have the same 
slope at A = 1 by design. 

We consider three individual cases for the anisotropic strains of (5.27) where the 
initial strong direction of the material a 0 aligns with each one of the three coordinate 
axes in Fig. 6.1. In the first case, 

a 0 = {l 0 0} T , (6.18) 


leading to 


\m-K 2 ) o 

o' 

and A ^ = 

riw-Ar 4 ) 

0 

o' 

0 0 

0 

0 

0 

0 

0 0 

0 


0 

0 

0 


which are associated with the anisotropic invariants (5.26) 

-^ = |(A 2 4- A x 2 ) an d h — ^(Ai + A! 4 ). 


(6.19) 


(6.20) 


In the second case, 

ao = {0 1 0} T , 

with 



O 

o 

o 


'0 0 

0" 

A (1) = 

0 |(a 2 -a 2 2 ) 0 

and A^ = 

0 |(a^-a 2 - 4 ) 

0 


0 0 0. 


0 0 

0 


( 6 . 21 ) 


( 6 . 22 ) 
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so that 


(6.23) 


h — \ (^2 + ^2 2 ) an d h — ^(^2 + ^ 2 4 )- 

And in the third case, 

a 0 = {0 0 1 } T , 

producing 


A™ = 


0 0 
0 0 
0 0 


0 

0 


m 


As' 2 ) 


and A (2) = 


0 

0 


| (A 3 - a 3 4 )J 


with 

h = | (A 3 + A 3 2 ) and I 5 = ^ (A 4 + A 3 4 ) . 


(6.24) 


(6.25) 


(6.26) 


In all three of these cases the non-zero anisotropic strain component equals the cor- 
responding isotropic strain component; they have the same strength. In fact, the 
coefficients appearing in the definitions of I 4 and 1 $ were chosen to specifically achieve 
this scaling. 

It is possible to run experiments where a 0 does not align with a coordinate axis, 
but care is required here because the ensueing deformations will not remain homoge- 
neous and they certainly will no longer be shear free. 


6.1.4 Stress Fields 

The actual state of stress is dictated by the boundary conditions imposed on the 
specimen. When tractions are present on a particular material plane they can lead 
to non-zero components of stress acting in that direction, otherwise these planes are 
stress (traction) free. 

For isotropic materials and transversely isotropic materials whose strong direction 
aligns with one of the three coordinate axes, the Cauchy stress, T, is characterized by 



'T 11 

0 

0 ‘ 


fi / 0 0 

T = 

0 

T 22 

0 

— 

0 /2/AlA3A 2 ,o 0 


0 

0 

T 33 . 


0 0 /3/AlA2A 3; o_ 


(6.27) 


which for incompressible materials (i.e., A 1 A 2 A 3 = 1 ) simplifies to 


T = 


Ai/i/A^o 

0 

0 


0 

Mf 2/^,0 

0 


0 

0 

Mfz/Azp 


(6.28) 


where /*, i = 1,2,3, are the contact forces at current time t > t 0 applied in the z th 
coordinate directions onto surfaces whose initial areas A i)0 (= fj,o4,od / j / k) are 
measured in the reference state of t 0 . 
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6.1.5 Special Cases 

Except for dilational compression, these experiments are considered to be done on 
(nearly) incompressible materials in the sense that their bulk moduli are order(s) in 
magnitude greater than their shear moduli. 

Simple tension 

Uniaxial tension/compression experiments done on isotropic materials have a defor- 
mation gradient and a state of stress with components that are described by 

AO 0 I \a = Xf/A 0 0 O' 

F — 0 A“V2 0 , T= 0 0 0, (6.29) 

0 0 A -1 /*] L 0 0 0 

where A (= Ai = £/l 0 ) is the applied stretch, with a being the resulting applied state 
of Cauchy stress. 

Biaxial tension 

When a sheet is stretched in orthogonal directions by amounts Ai and A 2 , the defor- 
mation that ensues is described by a deformation gradient and Cauchy stress whose 
matrix representations are 

Ai 0 0 <J\ — Ai/i/b4i,o 0 0 

F = 0 A 2 0 , T — 0 o-i — A 2 / 2 /A 2 ,o 0 , (6.30) 

_0 0 A^A^J [ ° 0 0_ 

where G\ and <r 2 are the corresponding applied stresses. Equibiaxial tension implies 
that Ai = A 2 . 

Pure shear 

When performed in a state of tension, this experiment is called pure shear; whereas, 
when performed in a state of compression, it is often called planar or channel com- 
pression. These experiments produce a state whose deformation gradient and stress 
fields have components given by 

'AO 0] for = A/Mo 0 O' 

F = 0 1 0 , T — 0 ? 0 , (6.31) 

0 0 A -1 J [ 0 0 0_ 

where it is under a condition of incompressibility that these kinematics conform with 
the classic notion of pure shear (in the sense of strains). Here A and A -1 are the 
applied and reaction stretches, respectively, while a and c are the resulting applied 
and reaction stresses, respectively. 

Dispite its name, pure shear is a shear-free deformation because the principal axes 
of stress and strain do not rotate, but remain orthogonal. Pure shear is therefore not 
a shearing deformation. Simple shear, defined in the next section, is a shearing 
deformation. 
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Dilational compression 

To be able to justify using the constraint of incompressibility, as we have done in 
the preceeding cases, one should experimentally demonstrate that the shear response 
of the material is orders of magnitude less stiff than its bulk response. A dilational 
compression experiment is usually performed to quantify bulk behavior. It is this 
experiment that earned Bridgman his Nobel Prize for Physics in 1946.* Here the 
deformation gradient and stress fields have matrix representations of 

'A 0 Ol [ct = f/A 0 0 O' 

F = 0 1 0 , T = 0 c 0 , (6.32) 

_o o ij [ o o c. 

where A is the applied stretch, a is the resulting applied stress, while c is the reaction 
stress, which Bridgman did not measure. 

This experiment does not, in general, impose a state of pure hydrostatic pressure. 
It is only when the bulk effects drown the shear effects that the resulting stress state 
becomes approximately isotropic (i.e., hydrostatic). 

6.2 Simple Shear 

This is a technologically important experiment for fluids and solids alike because the 
eigenvectors of stress and strain rotate in the body, yet the stress and strain fields 
themselves remain spatially uniform throughout the sample being tested. In this 
experiment, a homogeneous material element in the shape of a unit cube is deformed 
into a parallelepiped of unit height and with square shear planes by displacing the top 
shear plane a horizontal distance s, as illustrated in Fig. 6.3. This causes a reactionary 
shear stress, T 12 , to set up along this shear plane. Furthermore, second-order effects 
gradually appear in the normal directions as an artifact of these kinematics; effects 
first studied by Poynting [88]. 

The angle of shear, 7, and the magnitude of shear, s, are two common measures 
of shear deformation that relate to one another through the expression 

s = tan7 = 7 + |7 3 + ^7 5 + ^7 7 + --- , (6-33) 

with 7 being typically refered to as the engineering shear strain. 

6.2.1 Kinematics 

The deformation describing simple shear locates an arbitrary particle fp in test sam- 
ple ® with a set of spatial coordinates X — say, X = (Xi, X 2 , X 3 ) — in rectangular- 
Cartesian coordinate system C in a reference state at t 0 . Later, this particle occupies 

‘Bridgman’s Nobel Prize citation reads: 

“For the invention of an apparatus to produce extremely high pressures, and for the 
discoveries he made therewith in the field of high pressure physics.” 
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Figure 6.3: Simple shearing of a material element. 


a different location with coordinates x — say, x = (rti, x 2 , £3) — in the same coordinate 
system C , but in the current state of t. These coordinates relate to one another via 
the motion 

xi = Xi + sX 2 , x 2 = X 2 and x 3 = X 3 , (6.34) 

so that the deformation gradient tensors, F and £ -1 , have components of 


'1 

s 

O' 


*1 

-s O' 

0 

1 

0 

and F 1 = 

0 

1 

0 

0 

0 

1 


0 

0 

1 _ 


whose rotation matrix in turn has components 


cos # sin # 
— sin 9 cos 9 
0 0 


0 

0 , 
1 


wherein 9 = tan a ( s /2)=tan x (itan7). 


(6.35) 


(6.36) 


The left-stretch tensor has a matrix representation of 


with inverse 


V~ 


(1 4- sin 2 9)/ cos 9 sin# 0 
sin 9 cos 9 0 

0 0 1 


cos 9 — sin 9 0 

— sin 9 (1 -I- sin 2 9) / cos 9 0 
0 0 1 


while the right-stretch tensor has components of 


U 


cos 9 sin 9 0 

sin# (1 + sin 2 #)/cos# 0 
0 0 1 


(6.37) 


(6.38) 


(6.39) 
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with inverse 


u - 1 = 


(1 + sin 2 9)/ cos 9 
— sin 9 
0 


— sin 9 0 
cos 9 0 

0 1 


(6.40) 


The richness of tensor fields for describing the kinematics of deformation is fully 
realized in simple shear making it a more informative experiment than those with 
shear-free deformation histories. 

The velocity gradient tensor, L, has values of 


L = 


0 Ds 0 
0 0 0 , 
0 0 0 


(6.41) 


leading to rate-of-deformation, D, and vorticity, W_, tensors whose matrix representa- 
tions are 

r 0 IDs ol r 0 \Ds 0] 

(6.42) 



0 

1 Ds 

o' 


0 

1 Ds 

o' 

D = 

| Ds 

0 

0 

and W — 

-1 Ds 

0 

0 


0 

0 

0 


0 

0 

0 


It is the simplistic form of the vorticity tensor that makes simple shear such an 
important experiment for material characterization. 

, 

The lower-fractal rate-of-deformation tensor, D, has a representation of 


al 

D 


0 \D«s O' 

\D*s \(D«s 2 -2sD°s) 0 , 
0 0 0 


(6.43) 


while the upper-fractal rate-of-deformation tensor, D , has a representation of 


of 

D 


\{2sD*s-D°s 2 ) \D°s O' 
\D«s 0 0 , 

0 0 0 


(6.44) 


for 0 < a < 1. The component |(2sD“s - D°s 2 ) goes to zero in the limit as a goes 

ct[ oc\ 

to one (1) from below, which it must in order for D and D to be consistent with Q_ 
in this limit. 

Simple shear is isochoric because det(F) = 1, implying that this deformation is a 
volume preserving process, independent of constitutive assumption. 


6.2.2 Deformation Fields 


In the Eulerian frame, the contravariant-like, Finger, deformation tensor, and the 
covariant-like, Cauchy, deformation tensor, 5 -1 , have components of 


1 -|- s 2 

s 

O' 


' 1 

— s 

O' 

s 

1 

0 

and B 1 = 

— s 

1 + s 2 

0 

0 

0 

1 


_ 0 

0 

1 


(6.45) 
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In the Lagrangian frame, the covariant-like, Green, deformation tensor, Q, and its 
contravariant-like inverse, C _1 , have matrix representations of 



'1 

s 

O' 


1 + s 2 -s 

o' 

c = 

s 

1 + s 2 

0 

and G 1 — 

—s 1 

0 


_0 

0 

1 


O 

o 

1 


that evolve as 



' 0 

Ds 

O' 


2s Ds 

-Ds 

O' 

DC = 

Ds 

2s Ds 

0 

and DC- 1 = 

-Ds 

0 

0 


0 

0 

0 


0 

0 

0 


and whose Caputo derivatives are 


' 0 

D-s 

O' 


’D?s 2 

-D«s 

O' 

D°s 

D?s 2 

0 

and D^G- 1 = 

—D°s 

0 

0 

0 

0 

0 


0 

0 

0 


(6.46) 


(6.47) 


(6.48) 


for 0 < a < 1. 

Simple shear is the simplest experiment wherein the four deformation tensors g, 
B~ l , C and C _1 have matrix representations B, B _1 , C and C ^ 1 that are all different. 


6.2.3 Strain Fields 

Now that the various kinematic and deformation tensors and their rates have been 
quantified for simple shear, it is a straightforward process to establish any strain or 
strain-rate field of interest. Here we provide the five strain fields used in our material 
models. 

The three isotropic invariants of (5.13) are determined to be 

h = \ (3 + s 2 ) , I 2 = |(3 + 4s 2 + s 4 ) and J 3 = 2, (6.49) 


whose associated strains (5.15) have components of 


£(!) = I (B -B" 1 ) = 




0 

0 , 
0 


(6.50) 


E (2) = |(B 2 -B- 2 ) 


"|s 2 (2 + s 2 ) |s(2 + s 2 ) 0 

|s(2 + s 2 ) -|5 2 (2 + s 2 ) 0 

0 0 0 


(l + |s 2 )E (1) , (6.51) 


and 

e = i(|F|-|F|- 1 ) =0, (6.52) 

the latter of which is zero because the deformation is isochoric, and therefore 


6 — Qo, 


(6.53) 
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Strain Measures in Simple Shear 



Shear Displacement, s 

Figure 6.4: A comparison of E^ and for simple shearing. 

in accordance with the conservation of mass. The fact that E ^ is proportional to 
E^ within a scalar factor is unique to the deformation of simple shear. In other 
deformations, like shear-free flows, E ^ and are distinct measures of strain. 

Figure 6.4 presents a comparison between the shear components of E^ 2 \ E W and 
E ^ during a deformation of simple shear. Contrasting these plots with experimental 
data will provide the modeler with an indication as to which strain fields to include in 
a model for that material. Like shear-free deformations, E^ is a linear strain measure 
while is a nonlinear strain measure. 

Even though we do not presently know how to derive £C/ 2 ) f r0 m a free-energy 
expression, we still present E^ 2 ' 1 in Fig. 6.4 for comparison purposes. Of particular 
note is the fact that E^ 2 ^ — >■ 1 as s — >■ oo. If stress is to be proportional to strain, then 
the strain field is characteristic of fluid-like behavior, not solid-like behavior, 
thereby lending additional support for our selection of the integrity basis listed in 
(5.9a & 5.23a) for isotropic and transversely isotropic solids. If we were developing 
a theory for viscoelastic fluids instead of one for viscoelastic solids, then seeking an 
integrity basis that would lead to strains 0 < n < 1, would be an appropriate 
exercise. 

As with shear-free extensions, we investigate those instances where the strong 
material direction ao aligns with one of the three coordinate axes of Fig. 6.3. In the 
case where o 0 aligns with the direction of shearing, there 

ao = {l 0 0} T , (6.54) 

and this direction remains fixed throughout the deformation. This produces a longi- 
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tudinal shearing along the fibers, which themselves lie in the plane of shearing. Here 
the anisotropic strain tensors of (5.27) have components 


A™ = 


\ s 

-I * 2 


0 0 


and 


A< 2 > = (1 + §s 2 )A (1) = 


0 |s(2 + s 2 ) O' 

|s(2 + s 2 ) -M2 + * 2 ) 0 


whose associated invariants (5.26) are 

I 4 = \ (2 + s 2 ) and I 5 = iq (2 + 4s 2 4- s 4 ) , 


(6.55) 


(6.56) 


(6.57) 


indicating that an anisotropic effect is being picked up by this mode of shearing. 

In the second case, the strong direction is initially normal to the shearing plane 
and therefore aligns with the coordinate direction 


a 0 = {0 1 0} T . 


(6.58) 


Unlike the prior mode of shearing, here the material direction of strength rotates in 
the body during the deformation to a new direction of 

a 0 = { s/VT+s 1 l/S/TT^ 2 0} T (6.59) 


In the presence of this cross-axis shearing, the anisotropic strain tensors have com- 
ponents of 


A< 4 > = 


|s 2 \s O' 
is 0 0 
0 0 0 


and 


A< 2 > 


|s 2 (2 + s 2 ) |s(2 + s 2 ) O' 
|s(2 + s 2 ) 0 0 


0 


0 


0 


= (1 + P 


associated with invariants 

I 4 = |(2 + s 2 ) and I 5 = jg (2 + 4s 2 -I- s 4 ) . 


(6.60) 


(6.61) 


(6.62) 


Consequently, longitudinal and cross-axis shearing have the same, first-order, aniso- 
tropic response, differing only in their second-order normal responses. Unlike shear- 
free extensions, where the strength of the isotropic and anisotropic strains are the 
same, here the anisotropic shear component has half the strength of the isotropic shear 
component. This means that the simple shear experiment has half the sensitivity to 
material anisotropy that the simple extension experiment has. 
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The third case, called transverse shearing, places the strong direction in the shear 
plane lying normal to the shearing direction, so that 

a 0 = {0 0 1} T (6.63) 

Like longitudinal shearing, here the strong direction remains fixed throughout the 
deformation. Unlike the other two shearing modes, this mode of simple shear leads 
to anisotropic strain tensors that do not contribute to the overall strain field in that 

A (1) = A (2) = 0, (6.64) 

because the invariants are fixed, 

h = 72 and / 5 = Vs, (6.65) 

implying that transverse shearing occurs in the plane of isotropy. 

6.2.4 Stress Fields 

For isotropic materials and for transversely isotropic materials where the strong di- 
rection ao aligns with any one of the three coordinate directions, the most general 
state of stress that can arise from simple shearing produces components for Cauchy 
stress of [68, pp. 62-64] 

Tn Ti 2 0 fn /A x fxi/Ai 0 

T = T 21 T 22 0 = hx/Ax /22M2 0 , (6.66) 

0 0 T33_ _ 0 0 fw/As^ 

where fx 2 /A 2 = f 2 i/Ax because of symmetry in stress (i.e. , because T^ 2 = T 2 i ). No- 
tation f mn denotes a force acting in the m direction on a material surface of area A n 
whose unit normal is in the n direction in the rectangular-Cartesian coordinate system 
C displayed in Fig. 6.3. The third normal stress, T 33 , will generally be zero valued ex- 
cept for incompressible materials where T 33 will equal — p with p denoting a Lagrange 
multiplier that has been introduced to secure an incompressible deformation. 

Incompressible materials 

If a material is incompressible (i.e., its bulk response is orders of magnitude more 
stiff than its shear response), then the isotropic constraint of (3.20b) that enforces 
this condition will not permit all three normal components of stress to be uniquely 
determined. For this reason, simple shear experiments done on (apparently) incom- 
pressible materials result in at most three, independent, stress-like quantities that 
can be measured via experiment; they are: 

r := Ti 2 = T 21 , ’fq := Tn — T 22 , \l/ 2 := T 22 — T33, (6.67) 

where r is the shear stress, \Hq is the first normal-stress difference, and #2 is the 
second normal-stress difference. Taking the difference between any two normal-stress 
components subtracts out the unknown contribution arising from p, the Lagrange 
multiplier. Data from shear experiments done on viscoelastic fluids are usually re- 
ported in terms of these three stress measures. 
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Chapter 7 

Bulk Material Models 


In this report we are primarily concerned with biological and synthetic, polymeric 
elastomers whose bulk moduli are at least an order in magnitude more stiff than their 
shear moduli. For the purpose of materials characterization, it is therefore reasonable 
to assume that the bulk and shear responses are uncoupled. However, for the purpose 
of numeric computation in finite elements, it is advantageous to treat these materials 
as being compressible, even if only slightly so. 

In this chapter elastic and viscoelastic constitutive formulae that quantify pres- 
sure in terms of dilatation are presented and then used to solve the boundary-value 
problem of Bridgman’s experiment (6.32). There is much that the reader can learn 
about constructing fractional-order, viscoelastic, material models by first studying 
the simpler one-dimensional case of isotropic bulk behavior. 

7.1 Elastic Response 

In his early studies done shortly after the first world war, Bridgman [11] used the 
method of least squares to fit his isothermal experimental data to the empirical for- 
mula 

X dV ~ dV() . i 2 3 (v -\\ 

6 := — — = -ap + bp-cp, (7.1) 

dV o 

where he reported fitted parameters for a, b and c instead of presenting plots of the 
raw data. Even then his peak pressures were well in excess of 1 GPa. In his later 
works, Bridgman presented raw data in graphical form; fitted expressions were not 
given. In the above formula, <5(fP; f 0 , t ) represents a definition for dilatation used in the 
linear theory of elasticity, which is different from Hencky’s definition (3.15) given by 
k\(fP; t 0 , t) (= ln(dVydVo)), and which is different from our definition (5.15b) given by 
e(fP;f 0 ,£) (= ^{dV/dVo — dV 0 /dV)). Scalar p(*p, t) (= -f tr^) denotes hydrostatic 
pressure and dV represents the volume of mass element fp whose gauge volume is 
dVo- In most cases, Bridgman had only to fit his data with a quadratic polynomial 
in order to obtain a quality approximation to the data. In all cases, a linear fit was 
inadequate. 

Hencky [50] derived a governing equation for pressure from thermodynamics (viz., 
S fP = —QodW/dA), and then proceeded to propose two constitutive equations, the 
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simpler being 


fp = -nA, (7.2) 

wherein k (> 0) is the bulk modulus, p is the hydrostatic pressure, and A is his 
measure for dilatation. We point out that ^ = e A from the conservation of mass. 
Equation (7.2) is linear; it is Hencky’s definition for dilatation that is non-linear. 
This is a good e xam ple of what we are striving to achieve in this document: simple 
constitutive relations whose complexity (non-linearity) lies within the definitions of 
the fields themselves. 

Hencky’s constitutive equation has a strain-energy density of 


® = 2^( In V /det( 2» 1 '2 ) )' 


S 2^( ln ( det S)' 


2£>o 

K a 2 

= — A 2 
%Qo 


(7.3) 


which is quadratic in structure, like the strain-energy density from classic elasticity. 
Substituting (7.3) into (5.3)*, and taking its trace, reproduces Hencky’s constitutive 
equation of (7.2). 

Expanding Hencky’s constitutive formula (7.2) for pressure as a power series in 
the measure of dilatation used in infinitesimal strain analysis (i.e. , in terms of S) gives 


p/K = -8 + |5 2 - fS 3 + f| <5 4 - 
The inverted form of this series is 


137x5 
60 


J 5 + 0(<5 6 ). 


S = -p/n + §(p/k) 2 - §(p/k) 3 + 0((p/k) 4 ), 


(7.4a) 


(7.4b) 


which allows Hencky’s formula to be expressed in the format of Bridgman’s polynomial 
(7.1), where coefficients b and c can now be expressed in terms of constant a after an 
assignment of n := 1/a. A tabulation of such comparisons is provided in Table 7.1, 
where perfect fits have reported ratios of one. Over half of these comparison ratios 
are within a factor of two (i.e., between l /2 and 2) which, given the fact that these 
are coefficients to second- and third-order terms in p/n, is remarkable testimony as 
to the accuracy of Hencky’s simple formula for describing the bulk, elastic, material 
response. Not one of these ratios has a negative value; hence, the predicted and 
observed curvatures in these pressure- volume plots are in agreement. 

The values that Bridgman [11] measured for bulk moduli in 1923, and which are 
reported in Table 7.1, are about five to ten percent larger, on average, than their 
accepted values of today (see, e.g., http://www.webelements.com). Unavoidable 


•T his derivation makes use of the following properties of determinants: 

^ = ItIi _1 and Ito'-tI = Ito 1 ! Ill = 
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Metal 

Name 

Atomic 

Number 

Bulk Modulus 
k (GPa) 

2 nd order 

V( 3/(2k 2 )) 

3 rtl order 
c/(8/(3« 3 )) 

Li 

Lithium 

3 

12 

0.86 

— 

Na 

Sodium 

11 

6.5 

0.82 

16 

Mg 

Magnesium 

12 

34 

1.5 

— 

A1 

Aluminum 

13 

78 

1.3 

— 

K 

Potassium 

19 

2.9 

1.5 

1.7 

Ca 

Calcium 

20 

17 

0.92 

— 

Fe 

Iron 

26 

180 

4.1 

— 

Co 

Cobalt 

27 

200 

4.8 

— 

Ni 

Nickel 

28 

190 

5.1 

— 

Cu 

Copper 

29 

150 

3.4 

— 

Ge 

Germanium 

32 

82 

2.4 

— 

Sr 

Strontium 

38 

12 

0.72 

— 

Mo 

Molybdenum 

42 

280 

6.6 

— 

Pd 

Palladium 

46 

200 

5.2 

— 

Ag 

Silver 

47 

110 

3.0 

— 

Cd 

Cadmium 

48 

52-72 

2.5 

— 

Sn 

Tin 

50 

60 

1.7 

— 

Sb 

Antimony 

51 

42 

1.9 

— 

Ce 

Cerium 

58 

29 

0.99 

— 

Ta 

Tantalum 

73 

210 

0.73 

— 

W 

Tungsten 

74 

330 

12 

— 

Pt 

Platinum 

78 

290 

9.3 

— 

Au 

Gold 

79 

180 

4.1 

— 

Pb 

Lead 

82 

43 

2.0 

— 

Bi 

Bismuth 

83 

36 

1.5 

— 

U 

Uranium 

92 

110 

1.8 

— 


Table 7.1: Bulk moduli, as reported by Bridgman [11], and ratios of Bridgman’s 
higher-order coefficients (7.1) to those predicted by Hencky’s formula (7.2). 
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frictional forces that were likely present in his testing apparatus, and which he could 
not measure, are the probable cause of these errors. This does not in any way belit- 
tle his experimental results, nor diminish the significance of his findings. After all, 
twenty-three years after publishing Ref. [11], Bridgman was awarded the Nobel prize 
in physics for the invention of his testing apparatus and for the discoveries he made 
therewith. 

We can conclude from the comparisons presented in Table 7.1 that Hencky’s for- 
mula (7.2) is a viable material model for representing bulk material behavior to ex- 
tremely large states of pressure. Because most engineering applications do not expe- 
rience such enormous volume changes, it is reasonable to consider an approximation 
of the above construction. 


7.1.1 Theory for pressure 


Taking the trace of (5.16a), which is our phenomenological theory for isotropic elas- 
ticity, produces the following constitutive equation for hydrostatic pressure, 

f P = X^(2R,itr| (1) +21J )2 tr|( 2 ) + 32IJ 3 e). (7.5a) 

For materials whose bulk moduli are many times larger than their shear moduli, like 
those we are interested in, this expression simplifies to 


given that W 3 — k/2q 0 , 


(7.5b) 


which is a quadratic approximation to Hencky’s constitutive equation (7.2). Figure 7.1 
demonstrates that |(det£ — det£ _1 ) ~ IndetF over the interval */ 2 < det£ < 3 / 2 . 
This is a huge range for volume change from a practical point of view, recalling 
that det£ = dV/dV 0 , and as such, the constitutive formulae of (7.5b) and (7.2) are 
essentially equivalent in the realm of present-day engineering applications. 

From the conservation of mass, the dilational strain measure of (7.5b) can be 
expressed as 


e = i(de.£-detr 1 ) = i(^-^) 

(dV - dV 0 ) (dV + dV 0 ) 


dV 0 2 dV 

dV - dV 0 

whenever dV ~ dV 0 , 


(7.6) 


dV Q 


so that for infinitesimal volume changes 

dV - dV 0 


f p = -k 


dV n 


or equivalently p = —k 


dV - dV 0 
dV ’ 


(7.7) 


which is in agreement with the classical theory of elasticity. 

There will be instances when non-elastic aspects of a pressure/dilatation response 
need consideration (e.g., in capillary flows of viscoelastic liquids [64]) but, for most 
applications, the assumption of an elastic bulk response will be adequate whenever 
compressiblity effects need to be addressed. 
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Dilational Strain Measures 



Figure 7.1: Comparison of strain measures used to represent dilatation in constitutive 
equations for hydrostatic pressure: Hencky’s theory (7.2), our theory (7.5b), and 
linear elasticity. 

7.2 Viscoelastic Response 

Because pressure does not relax to zero at a fixed dilation, even in infinite time, fluids 
and solids both behave like solids in their bulk response. In what follows, we seek 
a solid-like representation within the constructs of Boltzmann’s [9] linear theory of 
viscoelasticity. 

Rewriting Hencky’s strain-energy density (7.3) as 

W= £- (in y^det^' 1 - 1)) > ( 7 - 8 ) 

substituting it into (5.4), and then taking its trace, leads to 

f p = -K (&(t) 4(0, t) + ^ m(t - t') A(t', t ) dt'^j , (7.9a) 

where 0(t - t') and M(t - t ') := d<5(t - t')/dt' are the bulk relaxation and memory 
functions, respectively, which are material properties to be acquired from experimen- 
tal data, and where A is Hencky’s definition (3.15) for dilatation. By envoking the 
additive and anti-symmetric property A(a,c ) = A(a,b ) + A(b,c), V b e [o, c], which 
applies to this strain definition, the above formula can also be expressed as 

fp=-K ^0(0) 4(0, t) - £ Tt(t - t') 4(0, t') dt'^j , (7.9b) 
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that after an integration by parts becomes 

fp = -K^S{t)A(0,0 + ) + J^ 0(£ — t') DA(t') dt'^j . (7.9c) 

All three formulae in (7.9) are mathematically equivalent. The first two expressions 
require dilatation to be continuous in time (i.e., A E C°), while the third expression is 
more restrictive in that it requires dilatation to be both continuous and differentiable 
in time (viz., A E C 1 ). 

It is also possible to construct a viscoelastic theory where dilatation is given in 
terms of pressure, but in our end application (finite elements), displacements are 
assigned to which forces respond, thereby placing (7.9) in the desired format. 

The above constitutive theory is, in and of itself, too general. To have engineering 
utility, one needs to assign appropriate functional forms to 0 and 9JI (i.e., one needs 
to create a model). There are a variety of ways that this can be done. One approach is 
to convert a known differential equation into the format of (7.9), which can always be 
done through, for example, an application of Laplace transform techniques, provided 
that the differential equation is linear. As an outcome of this procedure, the bulk 
relaxation and memory functions can be determined for that particular differential 
equation. 

7.2.1 Voigt solid 

The simplest viscoelastic solid of the differential type that one can consider is the 
Voigt solid, 

^-p(t) = -k(1 + pD)A(t), (7.10a) 

satisfying a homogeneous initial condition. The two material constants in the model 
are k and p. A bar is placed over the viscoelastic material constant to designate that 
it is associated with bulk deformation. The Voigt solid has a relaxation function given 
by 

0(£-t') = (l + /5)$(t-0, (7.10b) 

where 5 is the Dirac delta distribution (function). Consequently, its memory function, 
9Jt(£ — £') = d 0(£ — t')/dt', is undefined. 

The Voigt solid is therefore dismissed on physical grounds. It predicts that sound 
waves travel with infinite speed. For this reason, if this constitutive equation is to be 
used, the constraint of a homogeneous initial condition must be adhered to. Beware, 
this material model can cause instabilities to arise that are artifacts of the model. 
These are not physical instabilities; rather, they are numeric instabilities. 

7.2.2 Kelvin solid 

The next simplest material model that can be used to describe bulk viscoelatic be- 
havior is the Kelvin (or standard viscoelastic) solid, 

(l + r£>)(f p)(t) = -k( 1 + pD)A(t), T$rPo+ = -k^A 0 +, (7.11a) 
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which can handle an inhomogeneous initial condition, and therefore predicts a finite 
speed for sound. The bulk relaxation function is determined to be 

<$(f - f') = 1 + exp(-(f-t')/f), 1 < 0(* - 0 < ?. (7- llb ) 

which has an associated memory function of 

Wl(t - t ') = ^ exp(— (£ - 0 < Wl(t — f) < (7.11c) 

Both material functions are monotonic and bounded on the positive real line whenever 
t > 0, t S [0,t], and p > r > 0. Here k (> 0) denotes the rubbery (quasi-static) 
bulk modulus; /t(p/r) (> k) represents the glassy (dynamic) bulk modulus; fi (> 0) 
is the characteristic, bulk, relaxation time; and p (> f) is the characteristic, bulk, 
retardation (or creep) time. 

The differential equation in (7.11a) and the integral equations of (7.9) with mate- 
rial functions (7.11b & 7.11c) are equivalent formulations. Whether one chooses the 
differential or integral form of a particular model for use in analysis is largely a matter 
of personal taste that may be swayed by factors like: the boundary-value problem 
being solved, and the solution (e.g., numerical) methods available for solving it. 

The limited data that are available for characterizing bulk viscoelasticity suggest 
that the bulk response of a material is not as viscoelastic as its shear response. For 
example, the bulk viscoelastic response of polyisobutylene [77] has a ratio for pj t 
(i.e., the ratio of dynamic to quasi-static moduli) that is less than 10; whereas, for 
the same material, the shear response has a ratio of dynamic to quasi-static moduli 
(e.g., p/r) that exceeds 10 4 [33]. 

For this material model, sound is predicted to travel at a speed of v s = \J up/ f g 
(with q denoting mass density). This is a direct consequence of the fact that there 
exists a finite, inhomogeneous, elastic, initial condition whose effective (i.e., dynamic) 
bulk modulus happens to be k(p/t). Measuring both v s and g is therefore a viable 
way for experimentally evaluating the dynamic bulk modulus, n(p/f), in isotropic 
compressible materials. The static bulk modulus, k , can be obtained from a dilational 
compression test, like what Bridgman used. The time constants p and f can then 
be extracted from their ratio p/r by performing, for example, an additional stress 
relaxation experiment in dilational compression using the same compression fixture 
used previously to establish k. 

Measuring wave speeds to quantify the dynamic bulk modulus works well for 
fluids, where there is typically one type of wave present. However, additional care is 
warrented when this approach is applied to solids, because solids propagate two types 
of waves: logitudinal (bulk) and transverse (shear). 

7.2.3 Fractional-order models 

It is a straightforward matter to extend (7.10 & 7.11) to fractional order. The greater 
flexibility afforded by fractional models over their classical counterparts, when it 
comes time to correlate with experimental data, is justification enough for us to 
consider such extensions. 
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Fractional Voigt solid 

The simplest fractional-order viscoelastic solid that one can consider is the fractional 
Voigt solid introduced by Caputo [12] which, for dilatation, takes on the form 

irp(t) = —k( 1 + p a D°)A(t), (7.12a) 


and satisfies a homogeneous initial condition. This model has an additional constant 
a that controls the rate of evolution. Equation (7.12a) has a relaxation function of 


= 1 + p d 


T(l-fi) (f-f') a ’ 


(7.12b) 


that can be extracted directly from the definition of Caputo’s derivative (1.8a). When 
differentiated, the memory function is found to be 


m(t - t') = p & 


a 1 

T(1 - a) it - f') 1+fi ‘ 


(7.12c) 


Material function © is weakly singular at the upper limit of integration where t' — t, 
while 27t is strongly singular. 

The fractional Voigt solid is also dismissed on physical grounds, because it too 
propagates a sound wave with infinite speed due to the weak singularity present in 
the relaxation function. 


Fractional Kelvin solid 

The simplest, fractional-order, differential equation that is physically admissible for 
representing a viscoelastic solid is the standard FOV solid of (2.2) introduced by 
Caputo and Mainardi [13] which, as it applies to bulk response, is given by 


(1 + f 5 £>“) (f p)(t) = —k(1 + p*D;)A(t), &p 0+ = -k (| Y zV, (7.13a) 

and it can satisfy an inhomogeneous initial condition. This model has a bulk relax- 
ation function of [75] 


®(*~0 = 1 + P -/ E a (-{(t-t')/f) a y 1 < 0(t - f) < f , (7.13b) 

obtained by taking the Laplace transform of (7.13a). Equation (7.13b) is a dimen- 
sionless version of (2.5), which is the relaxation modulus for an FOV solid. Associated 
with © is the memory function 9 Jl, which is computed to be t 




t-t' 


, 0 <m(t- t') < oo. (7.13c) 


^Confer with Mainardi and Gorenflo [76] or Podlubny [86, pp. 21-22], for example, for a listing 
of derivatives and other useful properties of the Mittag-Leffler function. 
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Like (7.11b & 7.11c), (7.13b & 7.13c) are both positive monotone-decreasing functions 
provided that t > o ,fe [o,t], p > t > 0, and 0 < a < 1; however, because 
— {(t — t')/f) 6 )/(t — t') x {t — t') Q ~ l as t' — > t, it is apparent that this memory- 
function becomes unbounded at time t exhibiting a weak singularity. Here E a (x) = 
E a ,i(x) is the Mittag-Leffler function in one parameter, a, which reduces to the 
exponential function e x whenever a = 1, while E a ^(x) is the Mittag-Leffler function 
in two parameters, a and (3, (cf. §1.5). Again, k (> 0) is the rubbery bulk modulus; 
now n(p /f) a (> k ) represents the glassy bulk modulus; parameters f and p are still 
the characteristic, bulk, relaxation and retardation times, respectively; and a is a new 
material constant accounting for the fractal order of evolution. Their values can be 
extracted from the same set of experiments used to quantify the constants of (7.11). 

Experimental data [103] suggests that the bulk time constants f and p are much 
smaller than their counterparts r and p for shear (bulk transients fade faster), and 
that the fractal order for bulk evolution a. is smaller than the fractal order for shear 
evolution a (bulk behavior is less rate sensitive), both being bound to the interval 
(0, 1]. In other words, the bulk response tends to be more solid-like, while the shear 
response is more fluid-like. 

Surface plots of y = E a> i(—x a ) and y — —E a ^(—x a )/x, which are representative 
of 0 and 9 It in (7.13), respectively, are presented in Figs. 7.2 & 7.3. In both plots 
0.001 < a < 1, while 0 < x in the first plot and 0.001 < x in the second plot. 
Figure 7.2 is characteristic of stress-relaxation curves, while Fig. 7.3 is indicative of 
the extent of remembrance of past events. These plots clearly illustrate the influence 
that the fractal order of evolution, a, has on the (material) response function, y. 

When a is very small, but still greater than zero, the normalized relaxation func- 
tion y = E a (—x a ) of Fig. 7.2 drops immediately from a value of y = 1 at x — 0 to 
a value of y & 1/2 at x = 0 + , from which it slowly asymptotes to zero as x goes to 
infinity. Strictly speaking, the Mittag-Leffler function is not defined for a — 0. This 
would be the elastic boundary where, if the Mittag-Leffer function were defined, it 
would have to exhibit a discontinuous jump so that at a — 0 the response would be 
y = 1 for all x > 0 (i.e., there can be no relaxation in the elastic limit). At the other 
boundary of a = 1, relaxation is smooth and occurs at an exponential rate (viz., 
Ei(— x 1 ) = e _I ). For all values of a that lie inbetween, the Mittag-Leffler function 
provides a smooth monotone transition from elastic to nearly classic viscoelastic. 

Examining the normalized memory function y — —E ai0 (—x a )/x of Fig. 7.3, when 
a is very small but still greater than zero, the response y is observed to behave like 
an impulse function indicating that the material has a perfect knownledge of the 
current state, but it has almost no recollection of even the most recent of past states. 
As a grows toward a value of unity, the memory function continues to maintain a 
nearly perfect knowledge of the current state (i.e., y is infinite at x — 0, but the 
strength of this singularity diminishes as a moves toward one), and to this, it adds 
some recollection of past events — albeit a memory that rapidly fades away with the 
passing of time. At the upper boundary of a — 1, memory is lost at an exponential 
rate which, cureously, is the slowest rate of memory loss over the interval 0 < a < 1. 
Because y — 00 at x = 0 for 0 < a < 1, this kernel (7.13c) representing the memory 
function 911 of the bulk constitutive equation (7.9) produces an integral that has a 
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Figure 7.2: A 3D surface plot of y = E a ^{-x a ), which appears in the relaxation 
function 0 of the standard FOV solid. 


-x 1 E a ,oH-x a L 



X 


3 


Figure 7.3: A 3D surface plot of y = -E afi {-x a )/x, which appears in the memory 
function Wl of the standard FOV solid. 
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weak singularity at the upper limit of integration, like the fractional derivative (cf. 
Eqn. 1.8a). 

Choosing the integral equation of (7.9) (with the bulk relaxation and memory func- 
tions being given by Eqns. 7.13b & 7.13c) over its equivalent, fractional, differential 
equation (7.13a) has an advantage of disguising the fractal nature of the constitutive 
equation. However, selecting 0 (7.13b) and 911 (7.13c) without a priori knowledge 
of the fractional differential equation (7.13a) would be unlikely, at least they have 
not appeared in the literature as such.* For this reason, our preferred method for 
constitutive construction is to begin with a fractional differential equation, ensuring 
that it possesses physically meaningful initial conditions, and to then convert it into 
an equivalent integral expression (i.e. , to solve the differential equation). 

Furthermore, using the integral formulation of (7.9a) with memory function 9J1 as 
the kernel of integraion, instead of employing its variant (7.9c) that uses the relaxation 
function 0 as its kernel, leads to a convolution integral whose recollection of past 
events fades much more rapidly, and as such, should be better suited for numerical 
analysis. 


Alternate model: The viscoelastic models presented in (7.9-7.13) are all based 

on Hencky’s definition for dilatation, A = lndet£, defined in (4.38). In contrast, 
our viscoelastic theory (5.20) introduces a dilational strain measure defined by e = 
|(det F-detF -1 ), as established in (5.15b), which is a second-order accurate approx- 
imate to Hencky dilatation. In fact, whenever the bulk modulus is much stiffer than 
the shear modulus for a given material, the trace of (5.20a) produces a hydrostatic 
pressure response of 

fp=-K^S 3 (t)e(0,t) + J^Tl 3 (t-t')e(t',t)dt , ^j when 22J 3 = (7.14a) 

wherein 

f , ^ 1 { detF detF*' 

^ = 2 ^det£> ~~ det£ 

with F = F(0,t) and F? = F(0,i'). Assigning to this constitutive equation the 
following material functions 

0 3 (f -H) = 1 + ft E ai (- ((f - t')/n ) a 3 ) , (7.14b) 



and 


m 3 {t - 1') = -ft 


£q 3 .o(-((f-OA-3r) 


(7.14c) 


*This is the first appearance of a memory function that is representative of a fractional-order 
viscoelastic model. Relaxation functions representative of fractional-order viscoelastic models have 
been in the literature since Caputo’s [12] and Mainardi’s [13] pioneering papers, but they have not 
appeared (e.g., phenomenologically) in the general viscoelastic literature, excluding that subset of 
papers which makes use of the fractional calculus. 
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results in an admissible material model for bulk viscoelasticity with constants a 3 , /? 3 
and t 3 . This is the bulk, viscoelastic, material model that we advocate using. 

It is easily verified that e(a,c) ^ e(a,b ) + e(b,c) for any b G (a, c); in words, 
dilatation e is not additive and anti-symmetric in its dependency upon state. This 
means that (7.14a) and 

f P=-k (& 3 (t) e(0, 0 + ) + £ 0 3 (t - t') dt'^j (7.15a) 

are not equivalent constitutive formulae; (7.14a & 7.15a) are distinct. However, the 
integral equation in (7.15a), employing the material functions of (7.14b &; 7.14c), and 
the fractional differential equation 

(X + T s “C«)(¥ P )(() = - K (l + ^*D,«)e(0 > (). jgrPo* =-k (^)°*e(0,0 + ), 

(7.15b) 

are equivalent formulae upon setting /3 3 = ((p 3 — r 3 )/r 3 ) a3 . Consequently, the integral 
form of our viscoelastic model for bulk response (7.14) and the fractional differential 
equation (7.15b) are different material models, which is why we introduced constant 
As in place of (but not equivalent to) p 3 . To what extent and under what conditions 
they do differ remains to be determined. 

7.3 Bridgman’s Experiment 

Dilational compression is described by the boundary- value problem of (6.32). The 
state of stress given in (6.32) becomes hydrostatic (i.e. , a — <; = f /A 0 , wherein / is 
the force applied over an area of Ao) whenever shear effects can be neglected (e.g., 
the bulk modulus is significantly larger than the shear modulus). 

Under these conditions, and when expressed in terms of the parameters measured 
in an actual experiment, the elastic constitutive equation (7.5b) becomes 

X T = ~ n ( 7 - 16 ) 

where stretch A ratios the current specimen length t to its gauge length l§. This is a 
material model in one constant — the bulk modulus k. 

Likewise, the viscoelastic constitutive equation (7.14a) is given by 

X T 0 = ~ K (® 3(t) s ( A ~ A_1 ) + 1 dt ' ) - ( 7 - 17 ) 

where A = £{t)/l o and \ t ' = £(t ')/£ 0 are the stretch ratios, and where the material 
functions 0 3 and 9Jt 3 are given by (7.14b & 7.14c), respectively. This model has four 
material constants: k, a 3 , A 3 and r 3 . 
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Appendix A 

Table of Caputo Derivatives 


For the convenience of the reader, we provide this appendix where we give some 
Caputo-type derivatives of certain important functions.* We do not strive for com- 
pleteness in any sense, but we do want to give at least the derivatives of the classical 
examples. 

Throughout this appendix, a will always denote the order of the Caputo-type 
differential operator under consideration. We shall only consider the case a > 0 and 
a N, where N := {1, 2, 3, . . .} while N 0 := {0, 1,2,.. .}. We use the ceiling function 
[a] to denote the smallest integer greater than (or equal to) a, and the floor function 
[aj (= [a] - 1) to denote the largest integer (strictly) less than a. Recall that for 
a G N, the Caputo differential operator coincides with the usual differential operator 
of integer order. 

Moreover, E Vyil denotes the Mittag-Leffler function of the two parameters /i and 
v > 0, given by (cf. [31, Chp. 18]) 

°°^ 

B -« = gl>FT7 0’ 


•0 is the Digamma function, given by 


^{x) = 


r'(x) 

r(*r 


1 F 1 and 2 F\ denote the usual hypergeometric functions, i.e., 


\F x (a\b]z) = 


r(ft) r' r(q + k) k 
mt^r(b + k)k\ ’ 


aGK, — 6 ^Nq, 


(sometimes called the Kummer confluent hypergeometric function [1, Chp. 13]), the 
power series being convergent for arbitrary complex z, and 


2 Fi(a,fe; c; z) 


r(c) r(o + fc)r(fr + fc) ^ 
r (°) r ( 6 )^ r {c + k)k\ ~ 


fl,i) 6 E, — N 0 , 


‘Tables of Riemann-Liouville integrals and derivatives can be found in various places in the 
literature (cf. e.g., Podlubny [86] or Samko et al. [102]). We do not repeat those results here. 
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(the Gauss hypergeometric function [1, Chp. 15]), in which case the power series 
converges for all complex z with \z\ < 1 and may be extended analytically into the 
entire complex plane with a branch cut along the positive real axis from +1 to +oo. 
(Note that in the formulae below, we only need to evaluate this function for negative 
values of z, so the branch cut for positive 2 gives no problems.) Finally, i = y/^1 is 
the imaginary unit. 

1. Let f(x ) = x 3 . Here we have to distinguish some cases: 

0 if j G N 0 and j < [a] , 

r« + 1 ) 


(£>;/)(x) = 1 


-x 3 a if j G No and j > [a] 


r {j + 1 - a)‘ 

or j £ N and j > |_aj ■ 

2. Let f(x) — (x + c ) 3 for arbitrary c > 0 and j G R Then 


(D?f)(x) 


r (j + 1) c?- m- 1 


T(j - LaJ)r(fal -a + 1) 

3. Let f(x) = x 3 lnx for some j > |_aj. Then 


a 2 F 1 ( 1, [a] - j; for] -a + 1; -x/c). 


(D*f){x) 


r (j-W) 


W / 

W M -kT(j - a + 1)* 

r 0 + 1 ) „k-c 


j-a 


+ 


-x k - |aj) - ip(j - a + 1) + \nx). 


T(j - a + 1)' 

4. Let f(x) = exp (jx) for some j G R Then 

(D?f)(x) = jW X W- a E lM _ a+1 (jx). 

5. Let f(x) = sin jx for some j G R Here again we have two cases: 

I J 2r( ( ra 1 ) , tt + 1) WliM-a + liW 

— iFi(l; [a] — a + 1; —ijx)] if [a] is even, 

ifa Wl=ral-a + l;v.) 

+ iFi(l; [a] —a + 1; —ijx)\ if [a] is odd. 

6. Finally we consider f(x) = cos jx with some j G R As in the previous example, 
we obtain two cases: 


(D a J)(x) = { 


,1“1 (_l)f a V2T.fal-a 


X 1 


iFi(l; [a] —a + 1; ijx) 


(D?f)(x) = { 


2T([a] - a + 1) 

+ iFi(l; [a] —a + 1; —ijx)] if ("a] is even, 

2r(W- a + D ^M-a + i;vx) 

— iFi(l; [a] — a + 1; —ijx)] if fa] is odd. 
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Appendix B 


Automatic Integration 


In §1.5.2 we had found the need to calculate two integrals numerically in order to 
complete the evaluation of the Mittag-Leffler function. In this appendix we outline a 
possible strategy for the solution of these integrals. Essentially we will follow the ideas 
explained in QUADPACK [84], The routines introduced in that book are the generally 
accepted standard when looking for efficient and reliable quadrature algorithms. They 
are in the public domain, but they can also be found in many commercially distributed 
software packages. The source code is written in FORTRAN77 and may be obtained, 
for example, from the URL http://www.netlib.org/quadpack. For integrands of 
the form encountered in §1.5.2, it turns out that the routine DQAG (Double precision 
Quadrature, Adaptive, General purpose) is the method of choice. Our description of 
this routine follows a top-down structure (i.e., we first explain the general strategy 
without giving much information on the details, and at a later stage we fill in those 
details). 


B.l The Fundamental Strategy 

The fundamental idea of automatic integration is the following. The user supplies the 
integrand function and the bounds for the interval of integration (i.e., the data that 
define the integral in question) and a desired accuracy (i.e., a bound for the relative 
or absolute error that he/she is willing to accept). The routine is then supposed to 
return either an approximation for the correct value of the integral that is sufficiently 
accurate to satify the user’s requirements or an error flag if it fails to find such an 
approximation. 

Typically the algorithms try to achieve this goal by sub-dividing the interval of 
integration in an adaptive way, thus concentrating more quadrature nodes in areas 
where the integrand is difficult to approximate. This leads to a structure indicated 
in Alg. B.l. 

In practice, one may also terminate the WHILE loop in this algorithm when; 
for example, too many sub-divisions have taken place thereby using up all available 
memory, or too much computing time has been consumed. In such cases the algorithm 
should return an error flag. 
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Algorithm B.l Automatic Integration 

GIVEN two real numbers a and b , a function / : [a, b] — > C, 
and some tolerance e > 0 THEN 

Calculate an approximation Q[f] for f(x)dx and an estimate e for the 
error f* f{x)dx - Q[f] 

Initialize a list L of the approximations obtained so far 
by L := {([a, b], Q[f], e)} 

WHILE error tolerance is not satified DO 

Take the interval from L with the largest error estimate 
and remove it from the list; 

Bisect this interval; 

Calculate approximations and error estimates for the two newly obtained 
subintervals; 

Add these subintervals, their corresponding approximations, 
and their error estimates to the list L 
END of while loop 

RETURN the sum of the approximations obtained for all the intervals. 


B.2 Approximation of the Integral 

In Aig. B.l we have left open some of the key details. Specifically, we have not 
said how the approximation of the integral itself will be performed, and we have not 
indicated how the required error estimates can be found. We now turn our attention 
to these two questions. 

The basic idea behind the solution is the concept of Gauss-Kronrod integration. 
That is, we calculate two different approximations for the integral, both of which are 
of the form 

n 

1 a j,nf \ x j,n)t 
3=1 

with suitably chosen values of n, a JjTl , and x hn (j — 1, ... ,n). In particular, we 
begin with a first approximation that is just a Gauss quadrature formula with rii 
nodes, which is the (uniquely determined) quadrature formula that gives the exact 
value of the integral whenever the integrand is a polynomial of degree not exceeding 
2 ni — 1. This formula has been thoroughly investigated. For a recent survey we 
refer to Ref. [10] and the references cited therein. Specifically, there is no quadrature 
formula with the same number of nodes that is exact for a larger class of polynomials. 
Moreover, as stated in [10], both theoretical and practical evidence suggest that this 
method is a very good one. 


B.3 Approximation of Error Estimates 

From rather general considerations, it is known that one, single, quadrature formula 
can only give an approximation, but not both an approximation and an error estimate. 
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Therefore, to derive also an error estimate, it is necessary to introduce a second 
quadrature formula. The heuristic argument is as follows. Assuming that the second 
quadrature formula is much more accurate than the first, then the difference between 
the two approximations will be a good approximation for the error of the cruder of 
the two approximations, and therefore, a rather conservative (and thus quite reliable) 
upper bound for the error of the finer of the two approximations. Hence, we need 
a second quadrature formula that is significantly better than the first (Gauss-type) 
formula. In view of the quality of the Gauss formula, this can only be achieved by 
selecting a formula with n 2 nodes, where n 2 > n i- In principle, one could use a Gauss 
formula again, but this would be very uneconomical because a Gauss formula with 
ni nodes would have at most one node in common with a Gauss formula with n 2 
nodes, and so almost none of the information gathered so far (i.e. , the function values 
f( x j,n i), whose calculation is typically the most computationally expensive part of the 
algorithm) could be reused. To overcome this difficulty, Kronrod [59, 60] suggested 
to construct a formula that is nowadays called the Kronrod extension of the Gauss 
formula or, shortly, the Gauss-Kronrod formula. His formula is based on the Gauss 
formula with ni nodes; it uses n 2 = 2ni + 1 nodes and is constructed according to 
the following criteria: 

• The ni nodes Xj ni of the Gauss method form a subset of the n 2 nodes Xj <n2 of 
the Gauss-Kronrod method. 

• The remaining nodes and the weights a 7) „ 2 of the of the Gauss-Kronrod scheme 
are determined in such a way that the resulting method is exact for all polyno- 
mials of degree not exceeding 3ni + 1. 

A recent survey on Gauss-Kronrod formulas is given in [29]. 

Algorithms for the concrete calculation of the required nodes and weights are 
available (cf. the survey papers mentioned above and the references cited therein). 
For our purposes, it is sufficient to use the tabulated values given in [84]. In particular, 
following the suggestions made there, we propose using the Gaussian method with 
15 points in combination with the 31-point Kronrod extension for computing the 
integral with function K mentioned in Alg. 1.4 (the monotonic part). In view of the 
nice smoothness properties of the integrand, this gives us an approximation with quite 
high accuracy without too much computational effort. For the other integral in that 
algorithm, with the oscillatory integrand function P, these formulae may be unable 
to follow the oscillations properly; thus, we propose using the 30-point Gauss method 
together with its 61-point Kronrod extension. This is essentially also the method 
suggested for oscillatory integrals in the NAG Library [81] (see the documentation 
for routine D01AKF). 
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Appendix C 

Table of Pade Approximates for 
Mittag-Leffler Function 


This appendix presents a scheme for fast computations of 


E a (-x a ), 0 < a < 1, x > 0, 


suitable for finite elements. This function appears in the relaxation functions of both 
the FOV fluid (2.1) and solid (2.2). 

The defining series (1.11) is used on the interval 0 < x < 0.1, and an asymptotic 
series (1.16) is used whenever x > 15, as in Alg. 1.4, but in contrast with Alg. 1.4, a 
Pade approximate (or rational polynomial) is applied inbetween where 0.1 < x < 15; 
specifically, 


E a (-x a ) 


E 


k - 0 T(l+afc) 

an+aiz+osx 2 

1+6ix+62X 2 +()3X— 0 3 

_ v 4 (~ x E ak 

2^k=i r(i-afc) 


0 < x < 0.1 
0.1 < x < 15 
x > 15 


(C.l) 


where the coefficients ao, «i, ^ 2 , , &2 an d i >3 are listed in Table C.l for values of a 

varying between 0.01 and 0.99 by increments of 0.01. Also listed are the r resultants 
from each nonlinear regression with the tolerance set sufficient to achieve accuracies 
of r 2 > 0.999, except for a > 0.9 where such tight r 2 values could not be achieved for 
the order of approximate employed. The Pade approximates were obtained from fits 
to exact values (within machine precision) for the Mittag-Leffler function (obtained 
from Alg. 1.4) at 150 evenly spaced grid points over the interval [0.1, 15]. The form 
of this rational polynomial was selected after a study of numerous tables for function 
approximation listed in the appendices of Hart et al. [48]. Our fits were constrained 
to be exact at the end points so that transitions will be smooth when crossing a 
boundary between solution domains. 

A number written as 0.123456(7) stands for 0.123456 • 10' = 1, 234, 560. 
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Table C.l: Coefficients for Pade approximates of Mittag-Leffler function. 


a 

a o 

ai 

0,2 

bi 

62 

i >3 

P 

0.01 

-.424129(3) 

-.695856(6) 

-.391944(6) 

-.138637(7) 

-.795778(6) 

-.207400(3) 

0.999371 

0.02 

-.291078(4) 

-.231705(7) 

-.121393(7) 

-.460100(7) 

-.250385(7) 

-.127923(4) 

0.999345 

0.03 

-.159594(4) 

-.883786(6) 

-.485177(6) 

-.174711(7) 

-.101562(7) 

-.832687(3) 

0.999381 

0.04 

-.401668(4) 

-.151329(7) 

-.678102(6) 

-.298958(7) 

-.144743(7) 

-.138260(4) 

0.999114 

0.05 

-.480619(4) 

-.145550(7) 

-.643881(6) 

-.286638(7) 

-.139751(7) 

-.171118(4) 

0.999129 

0.06 

-.358700(4) 

-.953746(6) 

-.452754(6) 

-.186815(7) 

-.996694(6) 

-.160873(4) 

0.999313 

0.07 

-.431271(4) 

-.958824(6) 

-.424475(6) 

-.187447(7) 

-.952144(6) 

-.175012(4) 

0.999215 

0.08 

-.104813(5) 

-.220562(7) 

-.109649(7) 

-.427837(7) 

-.248747(7) 

-.594850(4) 

0.999401 

0.09 

-.113182(5) 

-.204162(7) 

-.931135(6) 

-.395568(7) 

-.215524(7) 

-.560193(4) 

0.999340 

0.10 

-.100534(5) 

-.165546(7) 

-.753470(6) 

-.319371(7) 

-.177305(7) 

-.527749(4) 

0.999361 

0.11 

-.771890(4) 

-.116952(7) 

-.530686(6) 

-.224658(7) 

-.126954(7) 

-.429168(4) 

0.999382 

0.12 

-.209949(4) 

-.331485(6) 

-.178194(6) 

-.628014(6) 

-.428531(6) 

-.183967(4) 

0.999349 

0.13 

-.201966(5) 

-.264628(7) 

-.118361(7) 

-.503944(7) 

-.292906(7) 

-.123451(5) 

0.999408 

0.14 

-.144185(5) 

-.177380(7) 

-.787251(6) 

-.336281(7) 

-.198184(7) 

-.923408(4) 

0.999420 

0.15 

-.358699(4) 

-.414163(6) 

-.180624(6) 

-.782040(6) 

-.463020(6) 

-.235432(4) 

0.999423 

0.16 

-.367384(5) 

-.387982(7) 

-.157708(7) 

-.732069(7) 

-.414017(7) 

-.218683(5) 

0.999358 

0.17 

-.834348(4) 

-.877610(6) 

-.381391(6) 

-.163991(7) 

-.101120(7) 

-.619007(4) 

0.999453 

0.18 

-.665920(4) 

-.661019(6) 

-.279508(6) 

-.123090(7) 

-.755801(6) 

-.495376(4) 

0.999452 

0.19 

-.880313(4) 

-.839225(6) 

-.352670(6) 

-.155465(7) 

-.970538(6) 

-.689880(4) 

0.999463 

0.20 

-.516937(4) 

-.463982(6) 

-.186660(6) 

-.857478(6) 

-.525412(6) 

-.391939(4) 

0.999446 

0.21 

-.181654(5) 

-.159929(7) 

-.655973(6) 

-.293380(7) 

-.187316(7) 

-.153968(5) 

0.999480 

0.22 

-.287739(5) 

-.241617(7) 

-.962412(6) 

-.441685(7) 

-.280701(7) 

-.244057(5) 

0.999477 

0.23 

-.283969(5) 

-.233523(7) 

-.935683(6) 

-.423970(7) 

-.277398(7) 

-.261404(5) 

0.999497 

0.24 

-.218950(5) 

-.174274(7) 

-.687682(6) 

-.314803(7) 

-.207911(7) 

-.208550(5) 

0.999503 

0.25 

-.184865(5) 

-.142847(7) 

-.555414(6) 

-.256670(7) 

-.171282(7) 

-.182554(5) 

0.999510 

0.26 

-.154930(5) 

-.116572(7) 

-.448241(6) 

-.208289(7) 

-.140906(7) 

-.159998(5) 

0.999520 

0.27 

-.252656(5) 

-.187367(7) 

-.720787(6) 

-.332362(7) 

-.230555(7) 

-.280571(5) 

0.999534 

0.28 

-.348286(5) 

-.250157(7) 

-.938597(6) 

-.441764(7) 

-.306951(7) 

-.392527(5) 

0.999538 

0.29 

-.223178(5) 

-.158819(7) 

-.597404(6) 

-.278265(7) 

-.198723(7) 

-.272201(5) 

0.999553 

0.30 

-.262152(5) 

-.180630(7) 

-.659350(6) 

-.315184(7) 

-.224602(7) 

-.321389(5) 

0.999556 

0.31 

-.110316(5) 

-.741001(6) 

-.263517(6) 

-.128658(7) 

-.919175(6) 

-.137449(5) 

0.999558 

0.32 

-.344887(5) 

-.231540(7) 

-.831618(6) 

-.398374(7) 

-.294522(7) 

-.473795(5) 

0.999578 

0.33 

-.253440(5) 

-.165943(7) 

-.579558(6) 

-.284133(7) 

-.210336(7) 

-.352888(5) 

0.999583 

0.34 

-.320037(5) 

-.206983(7) 

-.712936(6) 

-.351975(7) 

-.264342(7) 

-.467058(5) 

0.999592 

0.35 

-.788648(5) 

-.506620(7) 

-.172863(7) 

-.854865(7) 

-.654223(7) 

-.122026(6) 

0.999604 

0.36 

-.289781(5) 

-.195379(7) 

-.707274(6) 

-.324061(7) 

-.268137(7) 

-.555978(5) 

0.999616 

0.37 

-.142051(5) 

-.971770(6) 

-.354821(6) 

-.159385(7) 

-.136569(7) 

-.302386(5) 

0.999617 

0.38 

-.816551(5) 

-.523578(7) 

-.175872(7) 

-.860767(7) 

-.705295(7) 

-.156056(6) 

0.999642 

0.39 

-.575871(5) 

-.348533(7) 

-.107617(7) 

-.574012(7) 

-.450994(7) 

-.991720(5) 

0.999625 

0.40 

-.173808(5) 

-.117277(7) 

-.408516(6) 

-.187857(7) 

-.167833(7) 

-.429558(5) 

0.999651 

0.41 

-.103477(5) 

-.686777(6) 

-.231234(6) 

-.109392(7) 

-.977034(6) 

-.258727(5) 

0.999700 

0.42 

-.195086(5) 

-.129617(7) 

-.429282(6) 

-.204698(7) 

-.185582(7) 

-.514821(5) 

0.999681 

0.43 

-.333165(5) 

-.221252(7) 

-.718492(6) 

-.346552(7) 

-.318194(7) 

-.922518(5) 

0.999693 

0.44 

-.473286(5) 

-.311742(7) 

-.981295(6) 

-.484992(7) 

-.447200(7) 

-.134275(6) 

0.999708 

0.45 

-.147598(5) 

-.982326(6) 

-.304968(6) 

-.151324(7) 

-.142189(7) 

-.447789(5) 

0.999719 

0.46 

-.152914(5) 

-.104712(7) 

-.326140(6) 

-.159246(7) 

-.154553(7) 

-.517302(5) 

0.999725 

0.47 

-.167037(5) 

-.114533(7) 

-.346912(6) 

-.172774(7) 

-.169118(7) 

-.587605(5) 

0.999740 

0.48 

-.142611(5) 

-.991368(6) 

-.293586(6) 

-.148061(7) 

-.147062(7) 

-.531760(5) 

0.999752 

0.49 

-.654706(5) 

-.455618(7) 

-.130849(7) 

-.675094(7) 

-.675481(7) 

-.253173(6) 

0.999766 

0.50 

-.188800(3) 

-.135370(5) 

-.383675(4) 

-.198223(5) 

-.202768(5) 

-.797517(3) 

0.999777 

0.51 

-.516448(5) 

-.373170(7) 

-.102716(7) 

-.541357(7) 

-.559798(7) 

-.228252(6) 

0.999790 

0.52 

-.392154(5) 

-.289687(7) 

-.780114(6) 

-.415703(7) 

-.437101(7) 

-.185908(6) 

0.999801 

0.53 

-.135449(5) 

-.106392(7) 

-.289899(6) 

-.150142(7) 

-.164334(7) 

-.749106(5) 

0.999802 

0.54 

-.953841(4) 

-.748349(6) 

-.193214(6) 

-.104887(7) 

-.114460(7) 

-.530462(5) 

0.999822 

0.55 

-.793386(3) 

-.639333(5) 

-.159875(5) 

-.886439(5) 

-.979799(5) 

-.469973(4) 

0.999835 

1 0.56 

-.599232(4) 

-.443865(6) 

-.952614(5) 

-.619108(6) 

-.646389(6) 

-.288690(5) 

0.999847 
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Table C.l: Coefficients for Pade approximates of Mittag-Leffler function. 


a 

do 

ai 

a 2 

bi 

&2 

6s 

r z 

0.57 

-.337393(5) 

-.286190(7) 

-.663678(6) 

-.388535(7) 

-.438767(7) 

-.223581(6) 

0.999859 

0.58 

-.463960(5) 

-.405168(7) 

-.899289(6) 

-.544254(7) 

-.620515(7) 

-.324201(6) 

0.999871 

0.59 

-.151791(4) 

-.125852(6) 

-.241507(5) 

-.169238(6) 

-.185486(6) 

-.904156(4) 

0.999886 

0.60 

-.502729(4) 

-.464569(6) 

-.919187(5) 

-.611566(6) 

-.705142(6) 

-.377940(5) 

0.999896 

0.61 

-.480103(5) 

-.461103(7) 

-.860143(6) 

-.600297(7) 

-.697998(7) 

-.378472(6) 

0.999907 

0.62 

-.671512(4) 

-.708363(6) 

-.131044(6) 

-.906252(6) 

-.108843(7) 

-.626832(5) 

0.999908 

0.63 

-.312780(5) 

-.328986(7) 

-.533551(6) 

-.418611(7) 

-.494518(7) 

-.268704(6) 

0.999928 

0.64 

-.567954(5) 

-.632239(7) 

-.945394(6) 

-.794766(7) 

-.947107(7) 

-.509518(6) 

0.999937 

0.65 

-.149230(4) 

-.182752(6) 

-.259499(5) 

-.226154(6) 

-.275373(6) 

-.151317(5) 

0.999939 

0.66 

-.433396(5) 

-.549988(7) 

-.664464(6) 

-.674519(7) 

-.815645(7) 

-.407908(6) 

0.999953 

0.67 

-.752321(4) 

-.109682(7) 

-.126030(6) 

-.132094(7) 

-.164296(7) 

-.843703(5) 

0.999948 

0.68 

-.503149(5) 

-.782828(7) 

-.716973(6) 

-.932883(7) 

-.115421(8) 

-.500161(6) 

0.999960 

0.69 

-.113475(5) 

-.204800(7) 

-.158167(6) 

-.240096(7) 

-.301897(7) 

-.117538(6) 

0.999955 

0.70 

-.277604(3) 

-.535397(5) 

-.226639(4) 

-.621910(5) 

-.768192(5) 

-.150446(4) 

0.999976 

0.71 

-.551694(5) 

-.158192(8) 

-.673086(6) 

-.178833(8) 

-.232888(8) 

-.532630(6) 

0.999919 

0.72 

-.310305(5) 

-.102892(8) 

0.663348(5) 

-.115056(8) 

-.146317(8) 

0.202771(6) 

0.99994 

0.73 

-.160766(4) 

-.449254(7) 

-.116371(6) 

-.482429(7) 

-.682061(7) 

-.111373(6) 

0.999719 

0.74 

0.341935(4) 

-.241906(7) 

-.669717(5) 

-.252295(7) 

-.376790(7) 

-.808538(5) 

0.999510 

0.75 

0.173763(4) 

-.103357(7) 

-.513330(6) 

-.106853(7) 

-.212490(7) 

-.787594(6) 

0.999166 

0.76 

-.198004(4) 

-.170209(6) 

-.546410(6) 

-.212224(6) 

-.727938(6) 

-.906762(6) 

0.999483 

0.77 

-.255574(4) 

-.399973(6) 

-.408010(6) 

-.461840(6) 

-.889107(6) 

-.736776(6) 

0.999574 

0.78 

-.107310(4) 

-.102865(6) 

-.104591(6) 

-.125079(6) 

-.207143(6) 

-.205284(6) 

0.999673 

0.79 

-.494450(4) 

-.463767(6) 

-.357502(6) 

-.563528(6) 

-.802633(6) 

-.767703(6) 

0.999707 

0.80 

-.218066(4) 

-.219047(6) 

-.128115(6) 

-.262855(6) 

-.339154(6) 

-.302753(6) 

0.999721 

0.81 

-.534378(4) 

-.539996(6) 

-.258539(6) 

-.645497(6) 

-.766988(6) 

-.675138(6) 

0.999723 

0.82 

-.675467(4) 

-.688012(6) 

-.275764(6) 

-.819312(6) 

-.908437(6) 

-.800006(6) 

0.999715 

0.83 

-.702848(3) 

-.690578(5) 

-.240789(5) 

-.824283(5) 

-.847142(5) 

-.779293(5) 

0.999699 

0.84 

-.609980(4) 

-.583751(6) 

-.176849(6) 

-.697471(6) 

-.667838(6) 

-.642815(6) 

0.999674 

0.85 

-.103053(5) 

-.915547(6) 

-.247199(6) 

-.110278(7) 

-.963517(6) 

-.101397(7) 

0.999639 

0.86 

-.774771(4) 

-.667011(6) 

-.156824(6) 

-.805409(6) 

-.651177(6) 

-.733497(6) 

0.999594 

0.87 

-.838904(4) 

-.693431(6) 

-.141730(6) 

-.840712(6) 

-.624128(6) 

-.763282(6) 

0.999538 

0.88 

-.447750(4) 

-.359964(6) 

-.629983(5) 

-.437447(6) 

-.299094(6) 

-.396021(6) 

0.999468 

0.89 

-.390787(4) 

-.296304(6) 

-.446099(5) 

-.362820(6) 

-.222076(6) 

-.331368(6) 

0.999382 

0.90 

-.724688(4) 

-.509096(6) 

-.654940(5) 

-.630388(6) 

-.335377(6) 

-.584487(6) 

0.999277 

0.91 

-.152262(5) 

-.100158(7) 

-.106340(6) 

-.125254(7) 

-.574996(6) 

-.117309(7) 

0.999150 

0.92 

-.150108(5) 

-.942028(6) 

-.792810(5) 

-.118669(7) 

-.467003(6) 

-.112366(7) 

0.998998 

0.93 

-.106495(5) 

-.626399(6) 

-.398859(5) 

-.797901(6) 

-.255226(6) 

-766267(6) 

0.998816 

0.94 

-.153906(5) 

-.847297(6) 

-371308(5) 

-.109242(7) 

-.268144(6) 

-.106481(7) 

0.998600 

0.95 

-.600332(4) 

-.311566(6) 

0.756057(4) 

-.406277(6) 

-.710353(5) 

-.401715(6) 

0.998343 

0.96 

-.445797(5) 

-.214492(7) 

-.119855(5) 

-.284119(7) 

-.270455(6) 

-.285641(7) 

0.998041 

0.97 

-.355171(5) 

-.159674(7) 

0.206708(5) 

-.214673(7) 

-.423690(5) 

-.219239(7) 

0.997685 

0.98 

-.145397(5) 

-.610793(6) 

0.190402(5) 

-.834188(6) 

0.464083(5) 

-.865369(6) 

0.997270 

0.99 

-.105996(5) 

-.416113(6) 

0.204340(5) 

-.577797(6) 

0.755603(5) 

-.608788(6) 

0.996786 
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